Displacement when clip raster by shapefile

Question

I run a script to clip raster by shapefile. Even with layers own the same crs (epsg: 4674) the output raster be displaced. I didn't find out in my script what is caused the error.

Here is my code:

"""Clip a raster image using a shapefile"""

import gdal, gdalnumeric
import shapefile
from PIL import Image, ImageDraw


def imageToArray(i):
    """
    Converts a Python Imaging Library array to a gdalnumeric image.
    """
    a=gdalnumeric.numpy.fromstring(i.tobytes(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a


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)


def cut_raster_array(raster, shp, output):
    # Load the source data as a gdalnumeric array
    srcArray = gdalnumeric.LoadFile(raster)

    # Also load as a gdal image to get geotransform (world file) info
    srcImage = gdal.Open(raster)
    geoTrans = srcImage.GetGeoTransform()

    # Use pyshp to open the shapefile
    r = shapefile.Reader("%s" % shp)

    # Convert the layer extent to image pixel coordinates
    minX, minY, maxX, maxY = r.bbox
    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)

    # Multi-band image?
    #Check this modification in script in: http://karthur.org/2015/clipping-rasters-in-python.html
    try:
        clip = srcArray[:, ulY:lrY, ulX:lrX]

    # Nope: Must be single-band
    except IndexError:
        clip = srcArray[ulY:lrY, ulX:lrX]


    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    # Map points to pixels for drawing the county boundary
    # on a blank 8-bit, black and white, mask image.
    pixels = []
    for p in r.shape(0).points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    # Create a blank image in PIL to draw the polygon.
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    # Convert the PIL image to a NumPy array
    mask = imageToArray(rasterPoly)

    # Clip the image using the mask
    clip = gdalnumeric.numpy.choose(mask, (clip, 0)).astype(gdalnumeric.numpy.float32)

    # Save clipping as tiff
    gdalnumeric.SaveArray(clip, "%s.tif" % output, format="GTiff", prototype=raster)
    return clip


# Raster image to clip
raster = "raster/sema/custos_it_sema_2015xmunicipio_uc_multiplicacao.tif"
# Polygon shapefile used to clip
shp = "vector/teste_corte/clip_numpyarray.shp"
# Name of clipped raster file(s)
output = "vector/teste_corte/teste_corte"

array_cropped = cut_raster_array(raster, shp, output)

Answer 1

I used to do this in a similar way, but have since found that using rasterio and geopandas to be a quite succinct and easy way to do this.

import rasterio
import rasterio.mask

sf = geopandas.open_file('/path/to/shapefile') #This needs to contain only a single geometry shape

with rasterio.open('/path/to/dataset') as src:
    out_image, out_transform = rasterio.mask.mask(src, sf.geometry, crop=True)
    out_meta = src.meta

#Out image is the shapefile cropped image, while out_transform and out_meta contain georeferencing information.

out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("masked_image.tif", "w", **out_meta) as dest:
    dest.write(out_image)

This broadly follows information in the rasterio docs and in this question I asked previously about how to do this with xarray's.

Answer 2

You seem to be going through a series of hoops to do this. ArcGIS already supports this function with the Extract by Mask tool. Just make a feature layer from your shapefile and run the tool. That should do what you need.