0

I'm trying to mask a raster with population data with my area of interest. I reprojected the original population raster according to the AOI's crs, but after masking and writing the new file, the crs is set back to the original's crs, even though both inputs have the right crs (which is 32644). Can anyone point me to my mistake? I know I could reproject again, but there has to be a simpler way, right?

raster_crs = CRS(AOI.crs).to_epsg()

open original raster

with rio.open(pop_path, mode='r+') as pop:
    transform, width, height = calculate_default_transform(
        pop.crs, raster_crs, pop.width, pop.height, *pop.bounds)
    kwargs = pop.meta.copy()
    kwargs.update({
        'crs': raster_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

CRS(pop.crs).to_epsg()

4326

reproject original raster to AOI's crs

    with rio.open(pop_proj_tif, 'w', **kwargs) as pop_proj:
        for i in range(1, pop.count + 1):
            reproject(
                source=rio.band(pop, i),
                destination=rio.band(pop_proj, i),
                pop_transform=pop.transform,
                pop_crs=pop.crs,
                pop_proj_transform=transform,
                pop_proj_crs=raster_crs,
                resampling=Resampling.nearest)

CRS(AOI.crs).to_epsg()

32644

CRS(pop_proj.crs).to_epsg()

32644

mask reprojected raster with AOI

with rio.open(pop_proj_tif) as pop_proj:
    pop_out_image, out_transform = rio.mask.mask(pop_proj, AOI.geometry, crop=True)
    pop_out_meta = pop.meta

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

writing masked raster gives wrong CRS

    with rio.open(pop_proj_clip_tif, "w", **pop_out_meta) as pop_dest:
        pop_dest.write(pop_out_image)

CRS(pop_dest.crs).to_epsg()

4326

5
  • Is the projected raster saved on disk before you open it at stage "mask reprojected raster with AOI"?
    – user30184
    Commented Aug 10, 2022 at 5:32
  • @user30184 that's what the with rio.open(pop_proj_tif, 'w', **kwargs) as pop_proj: statement does, then the reproject writes directly to the output band destination=rio.band(pop_proj, i).
    – user2856
    Commented Aug 10, 2022 at 5:36
  • Could you still check what is the CRS when you do with rio.open(pop_proj_tif) as pop_proj: at stage "mask reprojected raster with AOI"?
    – user30184
    Commented Aug 10, 2022 at 5:44
  • 1
    Please provide a actual runnable code sample, including imports (so we can just substitute our data). Your code fails immediately as CRS, AOI, rio are not defined.
    – user2856
    Commented Aug 10, 2022 at 9:01
  • @user2856 I was wondering how to deal with this and couldn't find an answer in the help section. The problem is that AOI itself is a long script, which I thought would be redundant. I suppose you'd just replace vector and raster with some samples...nevertheless I should have added the imports. Will do next time, sorry
    – Danny
    Commented Aug 11, 2022 at 3:47

1 Answer 1

1

You are copying pop.meta instead of pop_proj.meta.

Replace:

# mask reprojected raster with AOI
with rio.open(pop_proj_tif) as pop_proj:
    pop_out_image, out_transform = rio.mask.mask(pop_proj, AOI.geometry, crop=True)
    pop_out_meta = pop.meta #<== wrong

With:

# mask reprojected raster with AOI
with rio.open(pop_proj_tif) as pop_proj:
    pop_out_image, out_transform = rio.mask.mask(pop_proj, AOI.geometry, crop=True)
    pop_out_meta = pop_proj.meta #<== this works

Also the reproject function signature doesn't have the following arguments - pop_transform, pop_crs, pop_proj_transform or pop_proj_crs and it's only a lucky accident that you didn't get a TypeError: reproject() got an unexpected keyword argument 'pop_transform' as the reproject function accepts a **kwargs argument for any additional arguments that should be passed to both the image to image transformer GDALCreateGenImgProjTransformer2() (for example, MAX_GCP_ORDER=2) and the GDALWarpOptions (for example, INIT_DEST=NO_DATA).

Use the following instead:

        reproject(
            source=rio.band(pop, i),
            destination=rio.band(pop_proj, i),
            src_transform=pop.transform,
            src_crs=pop.crs,
            dst_transform=transform,
            dst_crs=raster_crs,
            resampling=Resampling.nearest)

Or just drop them completely as they will be read from the source and dest bands anyway:

        reproject(
            source=rio.band(pop, i),
            destination=rio.band(pop_proj, i),
            resampling=Resampling.nearest)
1
  • Perfect, thank you so much. That saves me a lot of time.
    – Danny
    Commented Aug 11, 2022 at 3:19

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