How to open and clip raster data more efficiently when your clip data are in a different CRS

[lwasser]

Hi there @snowman2
i am working through a workflow where we teach students how to open and clip landsat scenes.
In this case I have a landsat scene and a geopandas GDF that are in different CRS's

I would like to do this - where I open the path to a single landsat band, and crop the data from_disk to conserve memory.

See code below which I believe can be more efficient

band_path  = 'path/to/tif/here'

#  Open the  data  -  i'd  rather  clip here so i can  take advantage of  from_disk
#  I'm  not  clipping because  it's  not  accepting  my crs  object  as i  thought it might  -  user  error?
band = rxr.open_rasterio(band_path, 
                         masked=True)  #.rio.clip(crop_bound_box[0],   ##  ideally i'd  chain   this all together  for  efficiency
#                                                from_disk=True, 
#                                                crs=crop_bound.crs).squeeze()

# Reproject the GDF data for clippping - here you NEED the CRS  from  the  raster opened  above
#  Could I somehow skip opening the raster just to get that  CRS?
if crop_bound.crs != band.rio.crs:
    crop_bound = crop_bound.to_crs(band.rio.crs)

# Turn bounding box (extent) into a shapely object for  clipping
crop_bound_box = [box(*crop_bound.bounds.loc[0])]

# Clip the data but given i opened it above i'm using more memory  / resources
band_crop = band.rio.clip(crop_bound_box, all_touched=True, from_disk=True)

What is the most efficient approach here? Is there a way to access the CRS of a tif file without opening it?
OR is there a way to use the crs= parameter in the clip function to force clip to reproject my Geometry?
I `do not wish to reproject the raster

Many thanks for any ideas.

OH NOTE: i am getting an error when i provide the CRS to the clip function. but i think before I go down that rabbit hole i'll see if there is a better approach first.

[snowman2]

The crs kwarg represents the CRS input of the geometries and will handle reprojecting the geometries for you. rioxarray uses GDAL and geopandas uses pyproj to transform the points, so the speed/results may vary between the two methods.

Note:
rioxarray only re-projects the geometries if the CRS are not the same:

rioxarray/rioxarray/raster_array.py

Lines 700 to 707 in 2049384

	if self.crs != crs:
	if LooseVersion(rasterio.__version__) >= LooseVersion("1.2"):
	geometries = rasterio.warp.transform_geom(crs, self.crs, geometries)
	else:
	geometries = [
	rasterio.warp.transform_geom(crs, self.crs, geometry)
	for geometry in geometries
	]

geopandas also skips re-projecting if the CRS are the same:
https://github.com/geopandas/geopandas/blob/fdf1fd7390752f37cbb1085fda2f7cf1c5ae4ebc/geopandas/array.py#L791-L793

[lwasser]

ok... so here is a more reproducible example. to achieve this I did use a rasterio context manager which is allowing me to grab the crs without reading the data - this actually is an ok approach i think but i'm testing it now on H4 data. but i'm still wondering if crs does what i think it does in .clipwhich is reproject the geometries.

Thank you for your patience!!

import os

import geopandas as gpd
from shapely.geometry import box
import earthpy as et
import rioxarray as rxr
os.chdir(os.path.join(et.io.HOME,
                      'earth-analytics',
                      'data'))


et.data.get_data('cold-springs-fire')

fire_boundary_path = os.path.join("cold-springs-fire",
                                  "vector_layers",
                                  "fire-boundary-geomac",
                                  "co_cold_springs_20160711_2200_dd83.shp")
fire_boundary = gpd.read_file(fire_boundary_path)

band_path = 'earthpy-downloads/landsat-coldsprings-hw/LC080340322016072301T1-SC20200306230205/LC08_L1TP_034032_20160723_20180131_01_T1_sr_band3.tif'

#  fire_boundary is a geodataframe
#  Turn that into an extent box for clipping
crop_bound_box = gpd.GeoDataFrame([1],
                                  geometry=[box(*fire_boundary.bounds.loc[0])])


# This returns  TypeError: 'Polygon' object is not  subscritable
data = rxr.open_rasterio(band_path,
                         masked=True).rio.clip(crop_bound_box.geometry,
                                               crs=crop_bound_box.crs,
                                               all_touched=True,
                                               from_disk=True).squeeze()

error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-17-4547bcbc3f88> in <module>
     27 
     28 # This returns  TypeError: 'Polygon' object is not  scritable
---> 29 data = rxr.open_rasterio(band_path,
     30                          masked=True).rio.clip(crop_bound_box.geometry,
     31                                                crs=crop_bound_box.crs,

~/opt/miniconda3/envs/earth-analytics-python/lib/python3.8/site-packages/rioxarray/raster_array.py in clip(self, geometries, crs, all_touched, drop, invert, from_disk)
    701                 geometries = rasterio.warp.transform_geom(crs, self.crs, geometries)
    702             else:
--> 703                 geometries = [
    704                     rasterio.warp.transform_geom(crs, self.crs, geometry)
    705                     for geometry in geometries

~/opt/miniconda3/envs/earth-analytics-python/lib/python3.8/site-packages/rioxarray/raster_array.py in <listcomp>(.0)
    702             else:
    703                 geometries = [
--> 704                     rasterio.warp.transform_geom(crs, self.crs, geometry)
    705                     for geometry in geometries
    706                 ]

~/opt/miniconda3/envs/earth-analytics-python/lib/python3.8/site-packages/rasterio/env.py in wrapper(*args, **kwds)
    383         else:
    384             with Env.from_defaults():
--> 385                 return f(*args, **kwds)
    386     return wrapper
    387 

~/opt/miniconda3/envs/earth-analytics-python/lib/python3.8/site-packages/rasterio/env.py in wrapper(*args, **kwds)
    606                                 param, full_kwds[param], inequality, version, reason))
    607 
--> 608             return f(*args, **kwds)
    609 
    610         return wrapper

~/opt/miniconda3/envs/earth-analytics-python/lib/python3.8/site-packages/rasterio/warp.py in transform_geom(src_crs, dst_crs, geom, antimeridian_cutting, antimeridian_offset, precision)
    103     """
    104 
--> 105     return _transform_geom(
    106         src_crs,
    107         dst_crs,

rasterio/_warp.pyx in rasterio._warp._transform_geom()

rasterio/_features.pyx in rasterio._features.OGRGeomBuilder.build()

TypeError: 'Polygon' object is not subscriptable

this works - but i reprojected first!

import rasterio as rio

# reproject the data 
band_path = 'earthpy-downloads/landsat-coldsprings-hw/LC080340322016072301T1-SC20200306230205/LC08_L1TP_034032_20160723_20180131_01_T1_sr_band3.tif'
#  fire_boundary is a geodataframe
#  Turn that into an extent box for clipping
crop_bound_box = gpd.GeoDataFrame([1],
                                  geometry=[box(*fire_boundary.bounds.loc[0])],
                                  crs=fire_boundary.crs)

with rio.open(band_path) as src:
    new_crs = src.crs

crop_bound_box_proj = crop_bound_box.to_crs(new_crs)

# This works after  reprojecting
data = rxr.open_rasterio(band_path,
                         masked=True).rio.clip(crop_bound_box_proj.geometry,
                                               crs=crop_bound_box_proj.crs, # probably don't  need  this now?
                                               all_touched=True,
                                               from_disk=True).squeeze()

data

[snowman2]

What is the output of rioxarray.show_versions()?

[snowman2]

I ask because (ref):

clipped = xds.rio.clip(crop_bound_box_proj.geometry, crop_bound_box_proj.crs)
# Note: If you have rasterio < 1.2 you will need convert the geometries to dict-like objects if the projection
# of the geometries differ from the raster. For example:
#
# from shapely.geometry import mapping
# geometries = crop_bound_box_proj.geometry.apply(mapping)

[lwasser]

Of course - -

rioxarray (0.3.0) deps:
  rasterio: 1.1.6
    xarray: 0.16.2
      GDAL: 3.1.4

Other python deps:
     scipy: 1.5.3
    pyproj: 3.0.0.post1

System:
    python: 3.8.6 | packaged by conda-forge | (default, Jan 25 2021, 23:22:12)  [Clang 11.0.1 ]
executable: /Users/leahwasser/opt/miniconda3/envs/earth-analytics-python/bin/python
   machine: macOS-10.14.6-x86_64-i386-64bit

looks like my rasterio version is a bit dated! let me update that now. i wouldn't have thought to look at that - thank you!

[lwasser]

ok that was it. thank you @snowman2 !
the code below now runs!
Oddly i commented out my geodataframe because rasterio was telling me that the geometry was not valid which i found more here about. instead i provided the shapely box directly and with the update to rasterio 1.2 it runs well! thank you!! i will be sure to tell the students to update versions and will update our environments accordingly. many thanks !!!

band_path = 'earthpy-downloads/landsat-coldsprings-hw/LC080340322016072301T1-SC20200306230205/LC08_L1TP_034032_20160723_20180131_01_T1_sr_band3.tif'
#  fire_boundary is a geodataframe
#  Turn that into an extent box for clipping
# crop_bound_box = gpd.GeoDataFrame([1],
#                                   geometry=[box(*fire_boundary.bounds.loc[0])],
#                                   crs=fire_boundary.crs)


# This returns  TypeError: 'Polygon' object is not  scritable
data = rxr.open_rasterio(band_path,
                         masked=True).rio.clip([box(*fire_boundary.bounds.loc[0])],
                                               crs=crop_bound_box.crs,
                                               all_touched=True,
                                               from_disk=True).squeeze()