1

I am trying to clip a Tiff raster with a shapefile. I follow the examples provided in rioxarray (here)

The link to raster file (here) and shp file (here)

import xarray as xr
import rioxarray
import geopandas as gpd
from shapely.geometry import box, mapping

#read the shp file
geodf = gpd.read_file('shp_file.shp')

# Assign the crs 
prj = [l.strip() for l in open('shp_file.prj','r')][0]
geodf.crs = prj

# Select one of the polygons
tmp = geodf[geodf['OBJECTID']==1]

da = xr.open_rasterio('raster_file.tif')

file_clipped = da.rio.clip(tmp.geometry.apply(mapping),tmp.crs)
file_clipped.rio.to_raster('file_clipped.tif')


However the resulted raster file (file_clipped.tif) has a different pixel size than the original one. The pixel size for the original raster is (463.3127165m, 463.3127165m) and for the clipped is (830.1019504m, 506.502885m). I am not sure what I do wrong. Is there any way to fix the raster pixel size when clipping?

EDIT1: Updated the codes and also uploaded the shp and raster files.

EDIT2: clipping the TIFF file with the rasterio.mask (in cases where polygons are separated, refer to the answers).

def tif_clip(tif_file, shp_file):
    from rasterio.mask import mask
    
    def getFeatures(gdf):
        """Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
        import json
        return [json.loads(gdf.to_json())["features"][0]["geometry"]]

    coords = getFeatures(shp_file)

    out_img, out_transform = mask(tif_file, shapes=coords, crop=True)

    out_meta = tif_file.meta.copy()
    #epsg_code = int(data.crs.data["init"][5:])

    out_meta.update(
        {
            "driver": "GTiff",
            "height": out_img.shape[1],
            "width": out_img.shape[2],
            "transform": out_transform,
            #"crs": pycrs.parse.from_epsg_code(epsg_code).to_proj4(),
            "crs": tif_file.crs,
        }
    )
    out_file = 'clip_raster.tif'
    with rasterio.open(out_file, "w", **out_meta) as dest:
        dest.write(out_img)

tif_file = rasterio.open('raster_file.tif')
geodf = gpd.read_file('shp_file.shp')
tmp = geodf[geodf['OBJECTID']==1]
tif_clip(tif_file,tmp)

This approach keeps the pixel size unchanged and also crop the extend to the extend of the shapefile. However, it still would be nice to be able to do the same thing with rioxarray for NetCDF files.

4
  • The pixel size should remain the same. What is the original pixel size versus the clipped? Can you share the files?
    – snowman2
    May 16, 2020 at 2:59
  • Thanks @snowman2 , I updated my question now it is fully reproducible.
    – Ress
    Jun 15, 2020 at 19:42
  • Looks like the files are not public - are you able to make them publicly available?
    – snowman2
    Jun 15, 2020 at 20:55
  • Sorry, I made them public. Also, I just noticed when I set the drop option in rio.clip to False the pixel size of the clipped raster is similar to the original. However, this option also keeps the extent of the clipped file similar to the original unlike clipping in let say Arcmap. Thanks
    – Ress
    Jun 15, 2020 at 21:09

1 Answer 1

1

From what I can tell, it is because you have two disconnected polygons in a multipolygon for the clip operation. This seems to produce undesired behavior and is not really supported.

clip chape

You can use the convex hull to get the boundary area of the shapes: enter image description here

If you use that to clip, you will maintain the resolution of the original raster in the clipped raster.

Suggested usage:

file_clipped = (
    da
    .rio.clip(tmp.geometry.convex_hull.apply(mapping), tmp.crs, drop=True)
    .rio.clip(tmp.geometry.apply(mapping), tmp.crs, drop=False)
)

EDIT: Fix has been added to rioxarray so that it now works without the workaround in 0.0.29+.

7
  • Yes, the problem is for cases where there are separate polygons. I updated my question with a function that uses rasterio.mask for clipping tiffs and it works fine. Is there any plan to add this rasterio capability to rioxarray? I mean its good to have it for NetCDF files through rioxarray. Using convex hull doesn't really work for me beacuse it alters the shapefile which impacts subsequent analyses. Thanks
    – Ress
    Jun 16, 2020 at 16:45
  • You should be able to only use the copy of the convex hull of the geometries. If you do that, it shouldn't impact the subsequent analysis.
    – snowman2
    Jun 16, 2020 at 16:57
  • Probably less of an issue when this is implemented: github.com/corteva/rioxarray/issues/115
    – snowman2
    Jun 16, 2020 at 16:58
  • Thanks, by the subsequent analysis I meant analysis on the clipped raster file. Clipping using convex hull make the area bigger than it actually is.
    – Ress
    Jun 16, 2020 at 18:07
  • I see what you mean. You could 1. Clip with convex hull drop=True and then 2. Clip with raw geometry drop=False. But, definitely adding complexity there.
    – snowman2
    Jun 16, 2020 at 19:24

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.