-
Notifications
You must be signed in to change notification settings - Fork 23
Open
Description
It is often not clear to me why the computation of a GeoBox fails. In this case, arr.odc.geobox simply returns None, without any further explanation of what data is missing. From the code, I see that the GeoState requires either gcp or transform to be set. Retrieving those is done through rather complicated functions. It would be great if the user could back-track a bit better where things went wrong and what information in the dataset is missing.
In my case, I assigned a coordinate system to a reprojected data array:
>>> print(arr)
<xarray.DataArray (eps: 1, ref_time: 1, y: 1000, x: 1000)> Size: 8MB
array([[[[0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[0. , 0. , 0. , ..., 0. ,
0. , 0. ],
...,
[0.00092668, 0.0008394 , 0.00075212, ..., 0. ,
0. , 0. ],
[0.00095773, 0.00087046, 0.00078318, ..., 0. ,
0. , 0. ],
[0.00098879, 0.00090151, 0.00081423, ..., 0. ,
0. , 0. ]]]])
Coordinates:
* eps (eps) int64 8B 0
* ref_time (ref_time) datetime64[ns] 8B 2025-11-28T06:00:00
lead_time timedelta64[ns] 8B 02:00:00
valid_time (ref_time) datetime64[ns] 8B 2025-11-28T08:00:00
lon (y, x) float64 8MB 5.95 5.955 5.959 5.964 ... 10.61 10.61 10.62
lat (y, x) float64 8MB 45.8 45.8 45.8 45.8 ... 47.81 47.81 47.81
spatial_ref int64 8B 0
Dimensions without coordinates: y, x
Attributes:
metadata: <earthkit.data.readers.grib.metadata.StandAloneGribMetadata...
parameter: {'centre': 'lssw', 'paramId': 500401, 'shortName': 'TOT_PRE...
geography: {'gridType': 'unstructured_grid', 'bitmapPresent': 0}
vref: geo
vcoord_type: surface
origin_z: 0.0
and the GeoState looks like this:
print(arr.odc._state)
GeoState(spatial_dims=('y', 'x'), crs_coord=<xarray.DataArray 'spatial_ref' ()> Size: 8B
array(0)
Coordinates:
lead_time timedelta64[ns] 8B 02:00:00
spatial_ref int64 8B 0
Attributes: (12/18)
crs_wkt: PROJCS["CH1903+ / LV95",GEOGCS["CH190...
semi_major_axis: 6377397.155
semi_minor_axis: 6356078.962818189
inverse_flattening: 299.1528128
reference_ellipsoid_name: Bessel 1841
longitude_of_prime_meridian: 0.0
... ...
longitude_of_projection_origin: 7.43958333333333
azimuth_of_central_line: 90.0
scale_factor_at_projection_origin: 1.0
false_easting: 2600000.0
false_northing: 1200000.0
spatial_ref: PROJCS["CH1903+ / LV95",GEOGCS["CH190..., transform=None, crs=CRS('PROJCS["CH1903+ / LV95",GEOGCS["CH1903+",DATUM["CH1903+",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6150"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4150"]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.9524055555556],PARAMETER["longitude_of_center",7.43958333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2056"]]'), geobox=None, gcp=None)
Still,
>>> arr.odc.geobox
None
Metadata
Metadata
Assignees
Labels
No labels