1

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)

enter image description here

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?

enter image description here

2

2 Answers 2

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. qgis save settings

3
  • Thanks for your answer. Unfortunately this doesn't solve the problem. The saved TIFF-file again is simply masked with opaque black. The masked out areas do not show as transparent.
    – petezurich
    Dec 5, 2022 at 20:18
  • If you want the image to have a transparency layer outside the GIS, you need to add an alpha channel to the image. I have modified my post to reflect this use case.
    – zeppeh
    Dec 7, 2022 at 13:18
  • Thanks for the additions. However, adding an alpha layer in QGIS doesn't seem to work as you describe. The saved file doesn't have an alpha. And using GIS software is not what I am looking for in this particular use case. But again – thanks for your effort. I appreciate that.
    – petezurich
    Dec 8, 2022 at 6:48
0
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)
6
  • Thanks for your answer! Unfortunately this doesn't work for me since I want to avoid writing a file to disk. In addition the resulting array out_img yields the same results as in my sample code. The areas that should be transparent aren't.
    – petezurich
    Dec 1, 2022 at 19:18
  • if you want to avoid writing a file to disk use GDAL.
    – Helios
    Dec 1, 2022 at 21:02
  • About your question, However, the cropped areas are black but not transparent. What do I do wrong? Use numpy where 0 == np.nan
    – Helios
    Dec 1, 2022 at 21:07
  • Thanks again for your thoughts. Using GDAL again means to create an intermediate file on disk. And I already tried np.where() with no success. This doesn't work with ìmshow()'s alpha parameter. It also creates new problems because a value of 0 can be within the bounds of the shape therefore in an area where I want the raster pixels unmasked.
    – petezurich
    Dec 2, 2022 at 5:43
  • because you need to create a copy of array and then filter the zeros like np.nan
    – Helios
    Dec 2, 2022 at 21:01

Your Answer

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

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