Crop raster with shapefile

Question

I want to crop a geo raster image with a geo shapefile using rasterio and geopandas. I struggle to properly create a mask so that the cropped areas actually display as transparent (and not black or any other color) when plotting with matplotlib.

How can I achieve this?

I am looking for a Python solution and I'd like to avoid creating an intermediate file on disk or use a GIS software. I looked through various SO posts and rasterio sample code to no avail. I tried manipulating numpy arrays, masked arrays and using a dataset MemFile, again with no success.

Sample code to illustrate the problem:

src = rasterio.open('test.tif')
shape = gpd.read_file("shape.shp")

fig, ax = plt.subplots(figsize=(6,6))
rasterio.plot.show(src, ax=ax)
shape.plot(color="red", ax=ax)

I would like to crop the raster so that pixels are only visible within the red star shape.

So I do this:

out_image, out_transform = rasterio.mask.mask(src, shape.geometry)

fig, ax = plt.subplots(figsize=(6,6))
rasterio.plot.show(out_image, transform=out_transform, ax=ax)

However, the cropped areas are black but not transparent. What do I do wrong?

Answer 1

Here is another approach to directly create a .tif mask with a transparent background using a combination of fiona and rasterio. However, this unfortunately writes the file to disk.

import fiona
import rasterio
from rasterio.mask import mask

def create_mask_from_shapefile(shapefile_filepath, corresponding_orthomosaic_filepath):

    # open shapefile
    with fiona.open(shapefile_filepath, 'r') as shapefile:
        shapes = [feature['geometry'] for feature in shapefile]

    # open rasterfile
    with rasterio.open(corresponding_orthomosaic_filepath, 'r') as src:
        out_image, out_transform = mask(src, shapes, crop=True) # setting all pixels outside of the feature zone to zero
        out_meta = src.meta

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

    output_file = 'mask.tif'

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

input_raster = 'orthophoto.tif'
input_shapefile = 'shape.shp'

create_mask_from_shapefile(input_shapefile, input_raster)

Here is one way to prevent the file from being written to disk via GDAL.

There is also an example in the GDAL/OGR cookbook for clipping a raster file with a corresponding shapefile. I personally have not used this yet, but have heard that it works.

Edit:
With qgis you can export a cropped TIFF file that has a transparent background outside of a GIS. I'm not sure yet how this can be achieved with gdal. Unfortunately, the file is also written to disk.

Answer 2

First If you have a shapefile change to geojson with this code:

import geopandas as gpd

esri_shapefile = gpd.read_file(r'input_esri_shapefile')

export_geojson = esri_shapefile.to_file(r'output_export_geojson', driver='GeoJSON')

Second open the raster and clip using rasterio and geojson

# Open the raster to clip

raster = rasterio.open(r'input_raster')

# Open geojson file

geojson = gpd.read_file(r'input_geojson')

# extract coordinates

def getFeatures(gdf):
    import json
    return [json.loads(gdf.to_json())['features'][0]['geometry']]

coords = getFeatures(geojson)
print(coords)

# Cut raster

out_img, out_transform = mask.mask(dataset=raster, shapes=coords, crop=True, nodata=0)
print(out_img.shape)

# Reshape 3D array to 2D array 

out_img = np.reshape(out_img, (out_img.shape[0]*out_img.shape[1], out_img.shape[2]))
print(out_img.shape)

# Copy metadata 

out_meta = raster.meta.copy()
print(type(out_meta))
print(out_meta)

# Update metadata

out_meta.update({'driver':'GTiff',
                 'width':out_img.shape[1],
                 'height':out_img.shape[0],
                 'count':1,
                 'dtype':'uint16', # float32, etc...
                 'crs':geojson.crs, 
                 'transform':out_transform,
                 'nodata':0})

# Save raster clip

with rasterio.open(fp=r'C:\Users\name\Desktop\results\clip_raster.tif', 
             mode='w',**out_meta) as dst:
             dst.write_band(1,out_img)