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.

Inspection showing displacement between polygon and raster

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.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
        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)

2 Answers 2


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:

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


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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.