1

I would like to first merge rasters and then clip the resulting raster using geopandas and rasterio. Here's what I came up with so far:

import rasterio
import glob
from rasterio import merge
from rasterio import mask
import geopandas

files = glob.glob("*.tif")
datasets = []

for file in files:
    ds = rasterio.open(file)
    meta = ds.meta.copy()
    datasets.append(ds)

out_data, out_transform = merge.merge(datasets)

meta["count"] = out_data.shape[0]
meta["height"] = out_data.shape[1]
meta["width"] = out_data.shape[2]
meta["transform"] = out_transform

for ds in datasets:
    ds.close()

with rasterio.open("merged.tif", "w", **meta) as ds:
    ds.write(out_data)

province = geopandas.read_file("province.shp")

with rasterio.open("merged.tif") as ds:
    out_data, out_transform = mask.mask(ds, shapes=province.geometry)

with rasterio.open("province.tif", "w", **meta) as ds:
    ds.write(out_data)

I opened up the result and it looks correct on QGIS.

What bugs me is that I have to write an intermediate file called merged.tif when I don't really need it in my workflow. All I care about is province.tif, yet rasterios heavy reliance on Python context managers and open file objects forces me to write this file.

I saw that rasterio also proposes a MemoryFile object which would in theory allow me to skip writing this file on my physical hard disk, but I'm nevertheless still required to use a python context manager, which just doesn't seem very pythonic and smells like poor design on the part of rasterio.

Is there any other alternative or workaround? Are there any higher level or better designed libraries to work with gridded geospatial data using python?

0

1 Answer 1

3

You can use rasterio.features.geometry_mask to mask your numpy array without writing a dataset.

import rasterio as rio
import glob
from rasterio import features
from rasterio import merge
import geopandas as gpd
import numpy as np

files = glob.glob("*.tif")
datasets = []

for file in files:
    ds = rio.open(file)
    datasets.append(ds)

    meta = ds.meta.copy()
    datasets.append(ds)

out_data, out_transform = merge.merge(datasets)

meta["count"] = out_data.shape[0]
meta["height"] = out_data.shape[1]
meta["width"] = out_data.shape[2]
meta["transform"] = out_transform

for ds in datasets:
    ds.close()

province = gpd.read_file("province.shp")

mask = features.geometry_mask(
            province.geometry,
            out_shape=out_data.shape[1:],
            transform=out_transform,
            all_touched=False,
            invert=True)

with rio.open("province.tif", "w", **meta) as ds:
    ds.write(np.ma.MaskedArray(out_data, mask))

enter image description here

enter image description here

enter image description here

2
  • @user32882 did this answer your question or is there something that I could improve in my answer?
    – user2856
    Commented Jan 9, 2023 at 8:15
  • @user32882 I never got around to actually testing this but I will assume your answer is correct.
    – user32882
    Commented Mar 25, 2023 at 17:12

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.