I have done this before with GDAL and Python, but was hoping to be able to use xarray to scale things up.

Let's say I have three datasets: 1. A large multiband geotiff where each band is a different date of the same parameter. Spatial resolution is r0 (say). ds_A say 2. A large multiband geotiff where each band is a value that is compared to a threshold and used to mask ds_A. Times and spatial extents are the same as ds_A, but resolution here is r1. Let's call this ds_B. 3. A vector layer with polygons. Let's call it vect. Lots of features in vect (>1e3, 1e4)

I want to (i) mask ds_A with ds_B and then extract (say) the spatial average of ds_A for all times for all the features in vect.

My approach so far is this, but this takes forever, in terms of memory and speed.

def do_mask(raster, vector, attr=None):
    """Either creates a 0/1 mask (if `attr` is `None`), or a 
    rasterises the label `vector` according to a field in `attr`, 
    0 otherwise."""
    g = gdal.Open(raster)
    w, h = g.RasterXSize, g.RasterYSize

    out_ds = gdal.GetDriverByName("MEM").Create("", w, h, 1, gdal.GDT_Byte)
    g = None
    if attr is None:
        gdal.Rasterize(out_ds, vector, burnValues=1, allTouched=True)
        gdal.Rasterize(out_ds, vector, attribute=attr, allTouched=True)
    return out_ds.ReadAsArray()

def mask_crop_stack(raster, vector, chunks, attr=None):
    """Create either a mask where a vector is True, or a set of labels
    from all features in a vector."""
    # Open the original dataset to get the coordinates and stuff
    # for rasterisation
    ds0 =  xr.open_rasterio(raster, chunks=chunks)
    # Create the mask or labels field.
    mask = do_mask(raster, vector, attr)
    # Create a DataArray with the mask/labels, and copy the metadata from the
    # original or target array.
    mask_ds = xr.DataArray(mask, coords={'x':ds.x, 'y':ds.y}, dims=list(ds.dims)[1:])
    # Throw away unmasked stuff to make the data more manageable
    cropped_ds = ds.where(mask_ds > 0, drop=True)
    return cropped_ds

# Open first file, and crop to vector area
ds_A = mask_crop_stack("ds_A.tif",
                       chunks={"band":1, "x":1024, "y":1024})
# Create labels DataArray
labels = mask_crop_stack("ds_A.tif",
                         chunks={"band":1, "x":1024, "y":1024},
# For some reason, the extents are slightly different
# So interpolate to match ds_A
labels = labels.interp(x=ds_A.x, y=ds_A.y, method="nearest")
# Open mask dataset, and crop spatially
ds_B = xr.open_rasterio("ds_B.tif",
                       "vect.geojson", chunks={"band":1, })
# Interpolate to match ds_A (change resolution from r1 to r0, say)
ds_B_rsmpl = ds_B.interp(x=ds_A.x, y=ds_A.y, method="nearest")
thresh = thresh_me
# Apply the mask in ds_B to ds_A
ds_A = ds_A.where(ds_B_rsmpl < thresh, other=np.nan)
# Apply a sanity check to ds_A, values need to be between 0 and 10000
ds_A = ds_A.where(ds_A > 0, other=np.nan)
ds_A = ds_A.where(ds_A < 10000, other=np.nan)

# Extract the time series and store in a dictionary
tseries = {label:ds_A.where(labels==label, other=np.nan).mean(axis=(1, 2)).to_pandas()
          for label in np.unique(labels) if label > 0}
# **Ideally**, I'd do processing pixel by pixel, 
# `apply_ufunc` might not work as processing outputs large number 
# of parameters per pixel.

The last line of the previous code takes forever, and often crashes the Python interpreter. I'm very new to xarray, so I'm guessing that I'm doing something wrong. My GeoTIFFs are tiled and blocky. I'm sure I'm missing something obvious here!


Your Answer

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

Browse other questions tagged or ask your own question.