Ytukyg

So I am working on a Arcpy toolbox in Arcgis that take DEM raster file into a certain treatement.
However I need to clip those image because the original ones are way too big and too long to process.
But thing is, Arcgis clipping tool changes the data type wich I cannot use then.

I started looking for codes to do that and it appears that GDAL librairy might help to clip a geotiff with a shapefile. Here is the code I followed with some minor changes to adapt to my 1-channel DEM: < https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html >

from osgeo import gdal, gdalnumeric, ogr, osr

from PIL import Image, ImageDraw

gdal.UseExceptions()





# This function will convert the rasterized clipper shapefile

# to a mask for use within GDAL.

def imageToArray(i):

    """

    Converts a Python Imaging Library array to a

    gdalnumeric image.

    """

    a=gdalnumeric.fromstring(i.tostring(),'b')

    a.shape=i.im.size[1], i.im.size[0]

    return a



def arrayToImage(a):

    """

    Converts a gdalnumeric array to a

    Python Imaging Library Image.

    """

    i=Image.fromstring('L',(a.shape[1],a.shape[0]),

            (a.astype('b')).tostring())

    return i



def world2Pixel(geoMatrix, x, y):

  """

  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate

  the pixel location of a geospatial coordinate

  """

  ulX = geoMatrix[0]

  ulY = geoMatrix[3]

  xDist = geoMatrix[1]

  yDist = geoMatrix[5]

  rtnX = geoMatrix[2]

  rtnY = geoMatrix[4]

  pixel = int((x - ulX) / xDist)

  line = int((ulY - y) / xDist)

  return (pixel, line)



#

#  EDIT: this is basically an overloaded

#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly

#  so we can pass these params off to CopyDatasetInfo

#

def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):

    ds = gdal.Open( gdalnumeric.GetArrayFilename(array) )



    if ds is not None and prototype_ds is not None:

        if type(prototype_ds).__name__ == 'str':

            prototype_ds = gdal.Open( prototype_ds )

        if prototype_ds is not None:

            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )

    return ds



def histogram(a, bins=range(0,256)):

  """

  Histogram function for multi-dimensional array.

  a = array

  bins = range of numbers to match

  """

  fa = a.flat

  n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)

  n = gdalnumeric.concatenate([n, [len(fa)]])

  hist = n[1:]-n[:-1]

  return hist



def stretch(a):

  """

  Performs a histogram stretch on a gdalnumeric array image.

  """

  hist = histogram(a)

  im = arrayToImage(a)

  lut = 

  for b in range(0, len(hist), 256):

    # step size

    step = reduce(operator.add, hist[b:b+256]) / 255

    # create equalization lookup table

    n = 0

    for i in range(256):

      lut.append(n / step)

      n = n + hist[i+b]

  im = im.point(lut)

  return imageToArray(im)



def main( shapefile_path, raster_path ):

    # Load the source data as a gdalnumeric array

    srcArray = gdalnumeric.LoadFile(raster_path)



    # Also load as a gdal image to get geotransform

    # (world file) info

    srcImage = gdal.Open(raster_path)

    geoTrans = srcImage.GetGeoTransform()



    # Create an OGR layer from a boundary shapefile

    shapef = ogr.Open("%s.shp" % shapefile_path)

    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )

    poly = lyr.GetNextFeature()



    # Convert the layer extent to image pixel coordinates

    minX, maxX, minY, maxY = lyr.GetExtent()

    ulX, ulY = world2Pixel(geoTrans, minX, maxY)

    lrX, lrY = world2Pixel(geoTrans, maxX, minY)



    # Calculate the pixel size of the new image

    pxWidth = int(lrX - ulX)

    pxHeight = int(lrY - ulY)



    clip = srcArray[ulY:lrY, ulX:lrX]



    #

    # EDIT: create pixel offset to pass to new image Projection info

    #

    xoffset =  ulX

    yoffset =  ulY

    print "Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset )



    # Create a new geomatrix for the image

    geoTrans = list(geoTrans)

    geoTrans[0] = minX

    geoTrans[3] = maxY



    # Map points to pixels for drawing the

    # boundary on a blank 8-bit,

    # black and white, mask image.

    points = 

    pixels = 

    geom = poly.GetGeometryRef()

    pts = geom.GetGeometryRef(0)

    for p in range(pts.GetPointCount()):

      points.append((pts.GetX(p), pts.GetY(p)))

    for p in points:

      pixels.append(world2Pixel(geoTrans, p[0], p[1]))

    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)

    rasterize = ImageDraw.Draw(rasterPoly)

    rasterize.polygon(pixels, 0)

    mask = imageToArray(rasterPoly)



    # Clip the image using the mask

    clip = gdalnumeric.choose(mask, 

        (clip, 0)).astype(gdalnumeric.uint8)





      clip[:,:] = stretch(clip[:,:])



    # Save new tiff

    #

    #  EDIT: instead of SaveArray, let's break all the

    #  SaveArray steps out more explicity so

    #  we can overwrite the offset of the destination

    #  raster

    #

    ### the old way using SaveArray

    #

    # gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)

    #

    ###

    #

    gtiffDriver = gdal.GetDriverByName( 'GTiff' )

    if gtiffDriver is None:

        raise ValueError("Can't find GeoTiff Driver")

    gtiffDriver.CreateCopy( "OUTPUT.tif",

        OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset )

    )



    # Save as an 8-bit jpeg for an easy, quick preview

    clip = clip.astype(gdalnumeric.uint8)

    gdalnumeric.SaveArray(clip, "OUTPUT.jpg", format="JPEG")



    gdal.ErrorReset()





if __name__ == '__main__':





    main( "shapefile", "DEM.tiff" )

However I got a "shape mismatch ValueError" :

<ipython-input-22-32e4e8197a02> in main(shapefile_path, raster_path, region_shapefile_path)

    166 

    167     # Clip the image using the mask

--> 168     clip = gdalnumeric.choose(mask,         (clip, 0)).astype(gdalnumeric.uint8)

    169 





/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in choose(a, choices, out, mode)

    399 

    400     """

--> 401     return _wrapfunc(a, 'choose', choices, out=out, mode=mode)

    402 

    403 



/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in _wrapfunc(obj, method, *args, **kwds)

     49 def _wrapfunc(obj, method, *args, **kwds):

     50     try:

---> 51         return getattr(obj, method)(*args, **kwds)

     52 

     53     # An AttributeError occurs if the object does not have



ValueError: shape mismatch: objects cannot be broadcast to a single shape

I tried to look around in the code where could it come from and I realized that this part is probably not working as it should:

minX, maxX, minY, maxY = lyr.GetExtent()

ulX, ulY = world2Pixel(geoTrans, minX, maxY)

lrX, lrY = world2Pixel(geoTrans, maxX, minY)

print("ulX, lrX, ulY, lrY : " , ulX, lrX, ulY, lrY) #image pixel coordinates of the shapefile

print(srcArray.shape)  #shape of the raster image 

clip = srcArray[ ulY:lrY, ulX:lrX] #extracting the shapefile zone from the raster image?

print(clip)

And returns:

('ulX, lrX, ulY, lrY : ', 35487, 37121, 3844, 5399)

(5041, 5041)

Seems that those indexes are at of bounds (but strangely python doesn't bother that much) and nothing is copied.
So I tried to change a bit the code to get the "real" pixel value corresponding to the area I wish to extract by using a shapefile corresponding to my total raster image:

#shapefile corresponding to the whole raster image

region_shapef = ogr.Open("%s.shp" % region_shapefile_path)

region_lyr = region_shapef.GetLayer( os.path.split( os.path.splitext( region_shapefile_path )[0] )[1] )



RminX, RmaxX, RminY, RmaxY = region_lyr.GetExtent()

RulX, RulY = world2Pixel(geoTrans, RminX, RmaxY)

RlrX, RlrY = world2Pixel(geoTrans, RmaxX, RminY)



#linear regression to find the equivalent pixel values of the clipping zone

pX = float(srcArray.shape[1])/(RlrX - RulX)

X0 = -(RulX*pX)

pY = float(srcArray.shape[0])/(RlrY - RulY)

Y0 = -(RulY*pY)



idXi = int(ulX*pX+X0)

idXf = int(lrX*pX+X0)

idYi = int(ulY*pY+Y0)

idYf = int(lrY*pY+Y0)



clip = srcArray[idYi:idYf, idXi:idXf]

print(clip)

Returns an array that really extracted values :

[[169.4 171.3 173.7 ... 735.6 732.8 729.7]

 [173.3 176.4 179.9 ... 734.3 731.5 728.7]

 [177.8 182.  186.5 ... 733.1 730.3 727.5]

 ...

 [ 73.3  77.5  83.  ... 577.4 584.9 598.1]

 [ 72.8  76.5  81.5 ... 583.1 593.  606.2]

 [ 71.3  74.7  79.  ... 588.9 599.1 612.3]]

Though I still have that goddamn :

<ipython-input-1-d7714555354e> in main(shapefile_path, raster_path, region_shapefile_path)

    170 

    171     # Clip the image using the mask

--> 172     clip = gdalnumeric.choose(mask,         (clip, 0)).astype(gdalnumeric.uint8)

    173 

    174     # This image has 3 bands so we stretch each one to make them



/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in choose(a, choices, out, mode)

    399 

    400     """

--> 401     return _wrapfunc(a, 'choose', choices, out=out, mode=mode)

    402 

    403 



/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in _wrapfunc(obj, method, *args, **kwds)

     49 def _wrapfunc(obj, method, *args, **kwds):

     50     try:

---> 51         return getattr(obj, method)(*args, **kwds)

     52 

     53     # An AttributeError occurs if the object does not have



ValueError: shape mismatch: objects cannot be broadcast to a single shape

Am I missing or misunderstanding something about it? I really start to lack of ideas so if someone got an idea please that would be appreciated.

Otherwise if you know another way to clip my DEM without altering it I'm fine as well.

You can achieve that in a much more simple way using gdal.Warp() and a shapefile as a cutline

from osgeo import gdal



input_raster = "path/to/yourDEM.tif" 

# or as an alternative if the input is already a gdal raster object you can use that gdal object

input_raster=gdal.Open("path/to/yourDEM.tif")

input_shape = "path/to/yourShapefile.shp" # or any other format

output_raster="path/to/outputDEM.tif" #your output raster file



ds = gdal.Warp(output_raster,

              input_raster,

              format = 'GTiff',

              cutlineDSName = input_shape, # or any other file format

              cutlineWhere="FIELD = 'whatever'" # optionally you can filter your cutline (shapefile) based on attribute values

              dstNodata = -9999) # select the no data value you like

ds=None #do other stuff with ds object, it is your cropped dataset. in this case we only close the dataset.