Skip to content

Commit

Permalink
Identified projection error as xee version
Browse files Browse the repository at this point in the history
  • Loading branch information
@kcartier-wri
kcartier-wri committed Sep 4, 2024
1 parent 6aab4c2 commit 7b2bbe6
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 53 deletions.
4 changes: 2 additions & 2 deletions .github/requirements.txt
View file
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
earthengine-api==0.1.408
geocube==0.4.2
geopandas==0.14.1
geopandas==0.14.4
rioxarray==0.15.0
odc-stac==0.3.8
pystac-client==0.7.5
pytest==7.4.3
xarray-spatial==0.3.7
xee==0.0.3
xee==0.0.15
utm==0.7.0
osmnx==1.9.3
dask[complete]==2023.11.0
Expand Down
24 changes: 2 additions & 22 deletions city_metrix/layers/layer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
import utm
import shapely.geometry as geometry
import pandas as pd
from pyproj import Transformer

MAX_TILE_SIZE = 0.5

Expand Down Expand Up @@ -274,18 +275,9 @@ def get_image_collection(
:param name: optional name to print while reporting progress
:return:
"""

# TODO may not be reprojecting!!!!

crs = get_utm_zone_epsg(bbox)

ds_orig = xr.open_dataset(
image_collection,
engine='ee',
scale=scale,
geometry=ee.Geometry.Rectangle(*bbox),
)

# Use a higher version for xee https://github.com/google/Xee/issues/118
ds = xr.open_dataset(
image_collection,
engine='ee',
Expand All @@ -299,18 +291,6 @@ def get_image_collection(
print(f"Extracting layer {name} from Google Earth Engine for bbox {bbox}:")
data = ds.compute()

crs_orig = ds_orig.attrs['crs']
crs_modified = ds.attrs['crs']
print('CRS: %s' % (crs_orig))
print('CRS: %s' % (crs_modified))

# x_val_orig = ds_orig.coords['lon'].values.min()
# x_val_modified = ds.coords['X'].values.min()
# print('Org_x: %s' % (x_val_orig))
# print('Modified_x: %s' % (x_val_modified))



# get in rioxarray format
data = data.squeeze("time").transpose("Y", "X").rename({'X': 'x', 'Y': 'y'})

Expand Down
51 changes: 22 additions & 29 deletions tests/test_layer_parameters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,98 +36,98 @@ def test_albedo_spatial_resolution():
class_instance = Albedo()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_alos_dsm_spatial_resolution():
class_instance = AlosDSM()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_average_net_building_height_spatial_resolution():
class_instance = AverageNetBuildingHeight()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_built_up_height_spatial_resolution():
class_instance = BuiltUpHeight()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_esa_world_cover_spatial_resolution():
class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_land_surface_temperature_spatial_resolution():
class_instance = LandSurfaceTemperature()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_nasa_dem_spatial_resolution():
class_instance = NasaDEM()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_natural_areas_spatial_resolution():
class_instance = NaturalAreas()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_ndvi_sentinel2_spatial_resolution():
class_instance = NdviSentinel2(year=2023)
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_tree_canopy_height_spatial_resolution():
class_instance = TreeCanopyHeight()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_tree_cover_spatial_resolution():
class_instance = TreeCover()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_urban_land_use_spatial_resolution():
class_instance = UrbanLandUse()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

def test_world_pop_spatial_resolution():
class_instance = WorldPop()
doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance)
assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution
assert 1==2
# assert 1==2

# def evaluate_resolution__property(obj):
# data = obj.get_data(BBOX)

def evaluate_resolution__property(obj):
is_valid, except_str = validate_layer_instance(obj)
if is_valid is False:
raise Exception(except_str)

# Double the default scale for testing
cls = get_class_from_instance(obj)

# Double the spatial_resolution for the specified Class
doubled_default_resolution = 2 * cls.spatial_resolution
obj.spatial_resolution=doubled_default_resolution

data = obj.get_data(BBOX)



expected_resolution = doubled_default_resolution
method, tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data)
print (method, tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max)
tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, y_min, y_max = estimate_spatial_resolution(data)
print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, y_min, y_max)

return expected_resolution, actual_estimated_resolution

Expand Down Expand Up @@ -170,34 +170,27 @@ def get_class_from_instance(obj):

def estimate_spatial_resolution(data):
y_cells = float(data['y'].size - 1)

x_min = None
y_min = data.coords['y'].values.min()
y_max = data.coords['y'].values.max()

try:
method = 'a'
crs_string = data.crs
if crs_string.startswith('PROJCS['):
crs = CRS.from_wkt(crs_string)
else:
crs = CRS.from_string(crs_string)
crs_unit = crs.axis_info[0].unit_name
crs_unit = crs.axis_info[0].unit_name
except:
method = 'b'
crs_string = data.rio.crs.data['init']
crs = CRS.from_string(crs_string)
crs_unit = crs.axis_info[0].unit_name

if crs_unit == 'metre':
diff_distance = y_max - y_min
elif crs_unit == 'degree' or crs_unit == 'm':
x_min = data.coords['x'].values.min()
diff_distance = get_distance_between_geocoordinates(y_min, x_min, y_max, x_min)
else:
raise Exception('Unhandled projection units: %s' % crs_unit)
raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs_string))

ry = round(diff_distance / y_cells)

tt = type(data)
return method, tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max
return tt, ry, crs, crs_unit, diff_distance, y_cells, y_min, y_max

0 comments on commit 7b2bbe6

Please sign in to comment.