Skip to content

Commit

Permalink
created tests and adjusted code to also update metadata if crs is dif…
Browse files Browse the repository at this point in the history 
…ferent to wgs84
  • Loading branch information
@AliceBalfanz
AliceBalfanz committed Aug 4, 2023
1 parent f461362 commit 4802593
Show file tree
Hide file tree
Showing 5 changed files with 173 additions and 44 deletions.
41 changes: 41 additions & 0 deletions test/core/gridmapping/test_base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
# noinspection PyProtectedMember
from xcube.core.gridmapping.helpers import _to_affine
from xcube.core.gridmapping.regular import RegularGridMapping
from xcube.core.new import new_cube

GEO_CRS = pyproj.crs.CRS(4326)
NOT_A_GEO_CRS = pyproj.crs.CRS(5243)
Expand Down Expand Up @@ -372,3 +373,43 @@ def test_ij_bboxes_from_xy_bboxes(self):
[360, 0, 720, 180],
[0, 340, 20, 360],
[-1, -1, -1, -1]], dtype=np.int64))

def test_to_dataset_attrs(self):
gm1 = TestGridMapping(**self.kwargs(xy_min=(20, 56),
size=(400, 200),
tile_size=(400, 200),
xy_res=(0.01, 0.01)))
transformed_gm = gm1.transform('EPSG:32633')

self.assertEqual(
{
'geospatial_bounds': 'POLYGON((20 56, 20 58.0, 24.0 58.0, 24.0 56, 20 56))',
'geospatial_bounds_crs': 'CRS84',
'geospatial_lat_max': 58.0,
'geospatial_lat_min': 56,
'geospatial_lat_resolution': 0.01,
'geospatial_lat_units': 'degrees_north',
'geospatial_lon_max': 24.0,
'geospatial_lon_min': 20,
'geospatial_lon_resolution': 0.01,
'geospatial_lon_units': 'degrees_east'
},
gm1.to_dataset_attrs())

self.assertEqual(
{
'geospatial_bounds': 'POLYGON((19.73859219445214 56.011953311099965, '
'19.73859219445214 57.96226854502516, 24.490858156706885 '
'57.96226854502516, 24.490858156706885 '
'56.011953311099965, 19.73859219445214 '
'56.011953311099965))',
'geospatial_bounds_crs': 'CRS84',
'geospatial_lat_max': 57.96226854502516,
'geospatial_lat_min': 56.011953311099965,
'geospatial_lat_resolution': 0.006391282582157487,
'geospatial_lat_units': 'degrees_north',
'geospatial_lon_max': 24.490858156706885,
'geospatial_lon_min': 19.73859219445214,
'geospatial_lon_resolution': 0.01443261852497102,
'geospatial_lon_units': 'degrees_east'},
transformed_gm.to_dataset_attrs())
138 changes: 121 additions & 17 deletions test/core/test_update.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray
import rioxarray # this is needed for adding crs to a dataset and used as rio

from test.sampledata import create_highroc_dataset
from xcube.core.update import update_dataset_attrs
Expand Down Expand Up @@ -185,35 +185,139 @@ def test_update_global_attributes_crs(self):
self.assertEqual('CF-1.7', ds2.attrs.get('Conventions'))
self.assertEqual('MIT', ds2.attrs.get('license'))
self.assertEqual('pipo', ds2.attrs.get('history'))
self.assertEqual(-20.0, ds2.attrs.get('geospatial_x_min'))
self.assertEqual(-18.0, ds2.attrs.get('geospatial_x_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_x_resolution'))
self.assertEqual('degrees_east', ds2.attrs.get('geospatial_x_units'))
self.assertEqual(12.0, ds2.attrs.get('geospatial_y_min'))
self.assertEqual(13.5, ds2.attrs.get('geospatial_y_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_y_resolution'))
self.assertEqual('degrees_north', ds2.attrs.get('geospatial_y_units'))
self.assertEqual(-20.0, ds2.attrs.get('geospatial_lon_min'))
self.assertEqual(-18.0, ds2.attrs.get('geospatial_lon_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_lon_resolution'))
self.assertEqual('degrees_east', ds2.attrs.get('geospatial_lon_units'))
self.assertEqual(12.0, ds2.attrs.get('geospatial_lat_min'))
self.assertEqual(13.5, ds2.attrs.get('geospatial_lat_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_lat_resolution'))
self.assertEqual('degrees_north', ds2.attrs.get('geospatial_lat_units'))
self.assertEqual('2018-06-01T00:00:00.000000000',
ds2.attrs.get('time_coverage_start'))
self.assertEqual('2018-06-06T00:00:00.000000000',
ds2.attrs.get('time_coverage_end'))
self.assertIn('date_modified', ds2.attrs)

ds1 = xr.Dataset(coords=coords_with_bnds)
ds1.rio.write_crs("epsg:4326", inplace=True,
grid_mapping_name="crs").reset_coords()
ds2 = update_dataset_attrs(ds1, global_attrs=output_metadata)

self.assertIsNot(ds2, ds1)
self.assertEqual('CF-1.7', ds2.attrs.get('Conventions'))
self.assertEqual('MIT', ds2.attrs.get('license'))
self.assertEqual('pipo', ds2.attrs.get('history'))
self.assertEqual(-20.0, ds2.attrs.get('geospatial_lon_min'))
self.assertEqual(-18.0, ds2.attrs.get('geospatial_lon_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_lon_resolution'))
self.assertEqual('degrees_east', ds2.attrs.get('geospatial_lon_units'))
self.assertEqual(12.0, ds2.attrs.get('geospatial_lat_min'))
self.assertEqual(13.5, ds2.attrs.get('geospatial_lat_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_lat_resolution'))
self.assertEqual('degrees_north', ds2.attrs.get('geospatial_lat_units'))
self.assertEqual('2018-06-01T00:00:00.000000000',
ds2.attrs.get('time_coverage_start'))
self.assertEqual('2018-06-05T23:00:59.000000000',
ds2.attrs.get('time_coverage_end'))
self.assertIn('date_modified', ds2.attrs)

def test_update_global_attributes_3031_crs(self):
num_x = 5
num_y = 5
num_times = 5

x_min = -2602050.
y_min = -1625550.

res = 100
res05 = res / 2
x = np.linspace(x_min + res05, x_min + num_x * res - res05, num_x)
y = np.linspace(y_min + res05, y_min + num_y * res - res05, num_y)
x_bnds = np.array([[v - res05, v + res05] for v in x])
y_bnds = np.array([[v - res05, v + res05] for v in y])
time = [pd.to_datetime(f'2018-06-0{i}T12:00:00') for i in
range(1, num_times + 1)]
time_bnds = [(pd.to_datetime(f'2018-06-0{i}T00:00:00'),
pd.to_datetime(f'2018-06-0{i}T23:00:59')) for i in
range(1, num_times + 1)]

coords = dict(time=(['time'], time),
y=(['y'], y),
x=(['x'], x))

coords_with_bnds = dict(time_bnds=(['time', 'bnds'], time_bnds),
y_bnds=(['y', 'bnds'], y_bnds),
x_bnds=(['x', 'bnds'], x_bnds),
**coords)

output_metadata = dict(history='pipo', license='MIT',
Conventions='CF-1.7')

expected_crs = "CRS84"

expected_bbox = (
'POLYGON((-121.99380296455976 -62.293519314442854, -121.99380296455976 '
'-62.299510171081884, -121.99083040333379 -62.299510171081884, '
'-121.99083040333379 -62.293519314442854, -121.99380296455976 '
'-62.293519314442854))')

ds1 = xr.Dataset(coords=coords)
ds1.rio.write_crs("epsg:3031", inplace=True,
grid_mapping_name="crs").reset_coords()
ds2 = update_dataset_attrs(ds1, global_attrs=output_metadata)

self.assertIsNot(ds2, ds1)
self.assertEqual('CF-1.7', ds2.attrs.get('Conventions'))
self.assertEqual('MIT', ds2.attrs.get('license'))
self.assertEqual('pipo', ds2.attrs.get('history'))
self.assertEqual(-121.99380296455976,
ds2.attrs.get('geospatial_lon_min'))
self.assertEqual(-121.99083040333379,
ds2.attrs.get('geospatial_lon_max'))
self.assertEqual(0.0005945389425789926,
ds2.attrs.get('geospatial_lon_resolution'))
self.assertEqual('degrees_east', ds2.attrs.get('geospatial_lon_units'))
self.assertEqual(-62.293519314442854,
ds2.attrs.get('geospatial_lat_min'))
self.assertEqual(-62.299510171081884,
ds2.attrs.get('geospatial_lat_max'))
self.assertEqual(0.0011981728729608676,
ds2.attrs.get('geospatial_lat_resolution'))
self.assertEqual(expected_crs, ds2.attrs.get('geospatial_bounds_crs'))
self.assertEqual(expected_bbox, ds2.attrs.get('geospatial_bounds'))
self.assertEqual('degrees_north', ds2.attrs.get('geospatial_lat_units'))
self.assertEqual('2018-06-01T00:00:00.000000000',
ds2.attrs.get('time_coverage_start'))
self.assertEqual('2018-06-06T00:00:00.000000000',
ds2.attrs.get('time_coverage_end'))
self.assertIn('date_modified', ds2.attrs)

ds1 = xr.Dataset(coords=coords_with_bnds)
ds1.rio.write_crs("epsg:3031", inplace=True,
grid_mapping_name="crs").reset_coords()
ds2 = update_dataset_attrs(ds1, global_attrs=output_metadata)

self.assertIsNot(ds2, ds1)
self.assertEqual('CF-1.7', ds2.attrs.get('Conventions'))
self.assertEqual('MIT', ds2.attrs.get('license'))
self.assertEqual('pipo', ds2.attrs.get('history'))
self.assertEqual(-20.0, ds2.attrs.get('geospatial_x_min'))
self.assertEqual(-18.0, ds2.attrs.get('geospatial_x_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_x_resolution'))
self.assertEqual('degrees_east', ds2.attrs.get('geospatial_x_units'))
self.assertEqual(12.0, ds2.attrs.get('geospatial_y_min'))
self.assertEqual(13.5, ds2.attrs.get('geospatial_y_max'))
self.assertEqual(0.25, ds2.attrs.get('geospatial_y_resolution'))
self.assertEqual('degrees_north', ds2.attrs.get('geospatial_y_units'))
self.assertEqual(-121.99380296455976,
ds2.attrs.get('geospatial_lon_min'))
self.assertEqual(-121.99083040333379,
ds2.attrs.get('geospatial_lon_max'))
self.assertEqual(0.0005945389425789926,
ds2.attrs.get('geospatial_lon_resolution'))
self.assertEqual('degrees_east', ds2.attrs.get('geospatial_lon_units'))
self.assertEqual(-62.293519314442854,
ds2.attrs.get('geospatial_lat_min'))
self.assertEqual(-62.299510171081884,
ds2.attrs.get('geospatial_lat_max'))
self.assertEqual(0.0011981728729608676,
ds2.attrs.get('geospatial_lat_resolution'))
self.assertEqual(expected_crs, ds2.attrs.get('geospatial_bounds_crs'))
self.assertEqual(expected_bbox, ds2.attrs.get('geospatial_bounds'))
self.assertEqual('degrees_north', ds2.attrs.get('geospatial_lat_units'))
self.assertEqual('2018-06-01T00:00:00.000000000',
ds2.attrs.get('time_coverage_start'))
self.assertEqual('2018-06-05T23:00:59.000000000',
Expand Down
4 changes: 1 addition & 3 deletions xcube/core/gen2/local/mdadjuster.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,4 @@ def _check_for_self_destruction(metadata: Dict[str, Any]):


def get_geospatial_attrs(gm: GridMapping) -> Dict[str, Any]:
# Note: check if force_geographic is wanted, here it is set for backwards
# compatibility
return dict(gm.to_attrs(force_geographic=True))
return dict(gm.to_dataset_attrs())
32 changes: 9 additions & 23 deletions xcube/core/gridmapping/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -612,17 +612,13 @@ def ij_bboxes_from_xy_bboxes(self,
ij_bboxes)
return ij_bboxes

def to_attrs(self,
force_geographic: bool = False) -> Mapping[str, Any]:
def to_dataset_attrs(self) -> Mapping[str, Any]:
"""
Get CF-compliant attributes of a dataset.
https://wiki.esipfed.org/Attribute_Convention_for_Data_Discovery_1-3#Recommended
Defined only for grid mappings with regular x,y coordinates.
:param force_geographic: Whether to force geographic CRS.
Defaults to false.
:return: dictionary with dataset coordinate attributes.
"""

Expand All @@ -649,24 +645,14 @@ def to_attrs(self,
lon_res = abs(lon_m2 - lon_m1)
lat_res = abs(lat_m2 - lat_m1)

if force_geographic:
geospatial_bounds_crs = 'CRS84'
geospatial_bounds = (f'POLYGON(('
f'{lon_min} {lat_min}, '
f'{lon_min} {lat_max}, '
f'{lon_max} {lat_max}, '
f'{lon_max} {lat_min}, '
f'{lon_min} {lat_min}'
f'))')
else:
geospatial_bounds_crs = self.crs.srs
geospatial_bounds = (f'POLYGON(('
f'{x1} {y1}, '
f'{x1} {y2}, '
f'{x2} {y2}, '
f'{x2} {y1}, '
f'{x1} {y1}'
f'))')
geospatial_bounds_crs = 'CRS84'
geospatial_bounds = (f'POLYGON(('
f'{lon_min} {lat_min}, '
f'{lon_min} {lat_max}, '
f'{lon_max} {lat_max}, '
f'{lon_max} {lat_min}, '
f'{lon_min} {lat_min}'
f'))')

return dict(
geospatial_lon_units='degrees_east',
Expand Down
2 changes: 1 addition & 1 deletion xcube/core/update.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ def update_dataset_spatial_attrs(dataset: xr.Dataset,
if update_existing or not all(coord_attr in is_in_attrs for coord_attr in
dataset.attrs):
# Update dataset with newly retrieved attributes
dataset.assign_attrs(gm.to_attrs())
dataset = dataset.assign_attrs(gm.to_dataset_attrs())
return dataset


Expand Down

0 comments on commit 4802593

Please sign in to comment.