Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
947 views
in Technique[技术] by (71.8m points)

geospatial - CRS error while clipping rioxarray to shapefile

I'm trying to clip a rioxarray dataset to a shapefile, but get the following error:

> data_clipped = data.rio.clip(shape.geometry.apply(mapping))
MissingCRS: CRS not found. Please set the CRS with 'set_crs()' or 'write_crs()'. Data variable: precip

This error seems straightforward, but I can't figure out which CRS needs to be set. Both the dataset and the shapefile have CRS values that rio can find:

> print(data.rio.crs)
EPSG:4326
> print(shape.crs)
epsg:4326

The dataarray within the dataset, called 'precip', does not have a CRS, but it also doesn't seem to respond to the set_crs() command:

> print(data.precip.rio.crs)
None
> data.precip.rio.set_crs(data.rio.crs)
> print(data.precip.rio.crs)
None

What am I missing here?


For reference, rioxarray set_crs() documentation - this shows set_crs() working on data arrays, unlike my experience with data.precip

My data, in case I have something unusual:

> print(data)
<xarray.Dataset>
Dimensions:      (x: 541, y: 411)
Coordinates:
  * y            (y) float64 75.0 74.9 74.8 74.7 74.6 ... 34.3 34.2 34.1 34.0
  * x            (x) float64 -12.0 -11.9 -11.8 -11.7 ... 41.7 41.8 41.9 42.0
    time         object 2020-01-01 00:00:00
    spatial_ref  int64 0
Data variables:
    precip       (y, x) float64 nan nan nan ... 1.388e-17 1.388e-17 1.388e-17
Attributes:
    Conventions:  CF-1.6
    history:      2021-01-05 01:36:52 GMT by grib_to_netcdf-2.16.0: /opt/ecmw...
> print(shape)
ID                      name                             orgn_name   geometry    
0                    Albania                             Shqip?ria   MULTIPOLYGON (((19.50115 40.96230, 19.50563 40...     
1                    Andorra                               Andorra   POLYGON ((1.43992 42.60649, 1.45041 42.60596, ...   
2                    Austria                            ?sterreich   POLYGON ((16.00000 48.77775, 16.00000 48.78252...     

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

This issue is resolved if the set_crs() is used in the same command as the clip operation:

data_clipped = data.precip.rio.set_crs('WGS84').rio.clip(shape.geometry.apply(mapping))

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...