0

I hope this doesn't get too confusing. I have two raster that overlap in a certain region. Like the following: enter image description here

What I would like to have, are two raster that only have values in the region where both rasters intersect, but also have values and no NA. Like the extent in the (really poorly drawn) image below:

enter image description here

What I did so far was using the gdaltindex command and compute a bounding box for each raster and intersect the, but than I get values for the red raster in the blue parts (and there I would like to have nodata in both resulting rasters)

enter image description here

I hope the explanation is sufficiently clear. Sorry about the terrible images... The solution could be on the command line with gdal, qgis, with rasterio in python, something in R....

Thanks a lot:)

3

I guess you can vectorize the part of your overlaping raster (the one on the right), then use this vector shape to crop the other raster (the red one).
Using rasterio library in Python, you can do the following:

import numpy as np
import rasterio
from rasterio import features
from rasterio.mask import mask

# the first one is your raster on the right
# and the second one your red raster
with rasterio.open('raster_with_no_data.tif') as src, \
        rasterio.open('raster_to_crop.tif') as src_to_crop:
    src_affine = src.meta.get("transform")

    # Read the first band of the "mask" raster
    band = src.read(1)
    # Use the same value on each pixel with data
    # in order to speedup the vectorization
    band[np.where(band!=src.nodata)] = 1

    geoms = []
    for geometry, raster_value in features.shapes(band, transform=src_affine):
        # get the shape of the part of the raster
        # not containing "nodata"
        if raster_value == 1:
            geoms.append(geometry)

    # crop the second raster using the
    # previously computed shapes
    out_img, out_transform = mask(
        dataset=src_to_crop,
        shapes=geoms,
        crop=True,
    )

    # save the result
    # (don't forget to set the appropriate metadata)
    with rasterio.open(
        'result.tif',
        'w',
        driver='GTiff',
        height=out_img.shape[1],
        width=out_img.shape[2],
        count=src.count,
        dtype=out_img.dtype,
        transform=out_transform,
    ) as dst:
        dst.write(out_img)
| improve this answer | |
2

Here's a method using R:

library(raster)
library(sf)
# Load rasters, make sure they are same resolution
r1 = raster("aridity.tif")
r2 = raster("landcover.tif")
r2_resamp = resample(r2, r1)

# Get raster extents and find intersection
r1_ext = st_as_sf(st_as_sfc(st_bbox(r1)))
r2_ext = st_as_sf(st_as_sfc(st_bbox(r2_resamp)))
r1_r2_intersect = st_intersection(r1_ext, r2_ext)

# Crop both rasters to the intersection
r1_crop = crop(r1, r1_r2_intersect, snap = "near")
r2_crop = crop(r2_resamp, r1_r2_intersect, snap="near")

# Make a mask where either raster is NA, and apply to each
na_mask = calc(stack(r1_crop, r2_crop), fun=sum, na.rm = FALSE)
r1_masked = mask(r1_crop, na_mask)
r2_masked = mask(r2_crop, na_mask)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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