From 386383dc70994987bb1b4e32644fc038c5215345 Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Tue, 24 Sep 2024 15:20:32 +0800 Subject: [PATCH 1/3] add GLAD LULC layers --- city_metrix/layers/__init__.py | 4 + city_metrix/layers/glad_lulc.py | 126 ++++++++++++++++++++++++++++++++ tests/test_layers.py | 20 +++++ 3 files changed, 150 insertions(+) create mode 100644 city_metrix/layers/glad_lulc.py diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py index 7e5e19e..a07449b 100644 --- a/city_metrix/layers/__init__.py +++ b/city_metrix/layers/__init__.py @@ -20,3 +20,7 @@ from .overture_buildings import OvertureBuildings from .nasa_dem import NasaDEM from .impervious_surface import ImperviousSurface +from.glad_lulc import LandCoverGlad +from.glad_lulc import LandCoverSimplifiedGlad +from.glad_lulc import LandCoverHabitatGlad +from.glad_lulc import LandCoverHabitatChangeGlad diff --git a/city_metrix/layers/glad_lulc.py b/city_metrix/layers/glad_lulc.py new file mode 100644 index 0000000..16d3234 --- /dev/null +++ b/city_metrix/layers/glad_lulc.py @@ -0,0 +1,126 @@ +import xarray as xr +import ee + +from .layer import Layer, get_utm_zone_epsg, get_image_collection + + +class LandCoverGlad(Layer): + """ + Attributes: + year: year used for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, year=2020, spatial_resolution=30, **kwargs): + super().__init__(**kwargs) + self.year = year + self.spatial_resolution = spatial_resolution + + def get_data(self, bbox): + data = get_image_collection( + ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}')), + bbox, + self.spatial_resolution, + "GLAD Land Cover" + ).b1 + + return data + + +class LandCoverSimplifiedGlad(Layer): + """ + Attributes: + year: year used for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, year=2020, spatial_resolution=30, **kwargs): + super().__init__(**kwargs) + self.year = year + self.spatial_resolution = spatial_resolution + + def get_data(self, bbox): + glad = LandCoverGlad(year=self.year, spatial_resolution=self.spatial_resolution).get_data(bbox) + # Copy the original data + data = glad.copy(deep=True) + + # Assign values based on ranges and conditions + data = data.where(glad != 0, 0) # If glad == 0, assign 0 + # If glad between 1 and 24, assign 1 + data = data.where((glad < 1) | (glad > 24), 1) + # If glad between 25 and 41, assign 2 + data = data.where((glad < 25) | (glad > 41), 2) + # If glad between 42 and 48, assign 3 + data = data.where((glad < 42) | (glad > 48), 3) + # If glad between 100 and 124, assign 4 + data = data.where((glad < 100) | (glad > 124), 4) + # If glad between 125 and 148, assign 5 + data = data.where((glad < 125) | (glad > 148), 5) + # If glad between 200 and 207, assign 6 + data = data.where((glad < 200) | (glad > 207), 6) + # data = data.where(glad != 254, 6) # Assign 6 if glad == 254 + data = data.where(glad != 241, 7) # If glad == 241, assign 7 + data = data.where(glad != 244, 8) # If glad == 244, assign 8 + data = data.where(glad != 250, 9) # If glad == 250, assign 9 + data = data.where(glad != 255, 10) # If glad == 255, assign 10 + + return data + + +class LandCoverHabitatGlad(Layer): + """ + Attributes: + year: year used for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, year=2020, spatial_resolution=30, **kwargs): + super().__init__(**kwargs) + self.year = year + self.spatial_resolution = spatial_resolution + + def get_data(self, bbox): + simplified_glad = LandCoverSimplifiedGlad(year=self.year, spatial_resolution=self.spatial_resolution).get_data(bbox) + # Copy the original data + data = simplified_glad.copy(deep=True) + + # Apply the conditions individually to avoid overwriting + data = data.where(simplified_glad != 0, 0) # Where glad == 0, set to 0 + # Where glad is between 1 and 6, set to 1 + data = data.where((simplified_glad < 1) | (simplified_glad > 6), 1) + data = data.where(simplified_glad < 7, 0) # Where glad >= 7, set to 0 + + return data + + +class LandCoverHabitatChangeGlad(Layer): + """ + Attributes: + start_year: baseline year for habitat change + end_year: target year for habitat change + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, start_year=2000, end_year=2020, spatial_resolution=30, **kwargs): + super().__init__(**kwargs) + self.start_year = start_year + self.end_year = end_year + self.spatial_resolution = spatial_resolution + + def get_data(self, bbox): + habitat_glad_start = LandCoverHabitatGlad(year=self.start_year, spatial_resolution=self.spatial_resolution).get_data(bbox) + habitat_glad_end = LandCoverHabitatGlad(year=self.end_year, spatial_resolution=self.spatial_resolution).get_data(bbox) + + # Class 01: Became habitat between start year and end year + class_01 = ((habitat_glad_start == 0) & (habitat_glad_end == 1)) + # Class 10: Became non-habitat between start year and end year + class_10 = ((habitat_glad_start == 1) & (habitat_glad_end == 0)) + + # Initialize the habitat change DataArray with zeros + habitatchange = xr.full_like(habitat_glad_start, 0) + # Set values 1 for class_01 (became habitat between start year and end year) + habitatchange = habitatchange.where(~class_01, other=1) + # Set values 10 for class_10 (became non-habitat between start year and end year) + habitatchange = habitatchange.where(~class_10, other=10) + + return habitatchange diff --git a/tests/test_layers.py b/tests/test_layers.py index d26a917..67ce921 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -9,6 +9,10 @@ EsaWorldCoverClass, HighLandSurfaceTemperature, ImperviousSurface, + LandCoverGlad, + LandCoverHabitatGlad, + LandCoverHabitatChangeGlad, + LandCoverSimplifiedGlad, LandsatCollection2, LandSurfaceTemperature, NasaDEM, @@ -60,6 +64,22 @@ def test_impervious_surface(): data = ImperviousSurface().get_data(BBOX) assert np.size(data) > 0 +def test_land_cover_glad(): + data = LandCoverGlad().get_data(BBOX) + assert np.size(data) > 0 + +def test_land_cover_habitat_glad(): + data = LandCoverHabitatGlad().get_data(BBOX) + assert np.size(data) > 0 + +def test_land_cover_habitat_change(): + data = LandCoverHabitatChangeGlad().get_data(BBOX) + assert np.size(data) > 0 + +def test_land_cover_simplified_glad(): + data = LandCoverSimplifiedGlad().get_data(BBOX) + assert np.size(data) > 0 + def test_land_surface_temperature(): data = LandSurfaceTemperature().get_data(BBOX) assert np.size(data) > 0 From 9df8d4b203ce1a7e4db236298d8b0b83325a66f3 Mon Sep 17 00:00:00 2001 From: Chris Rowe Date: Thu, 31 Oct 2024 18:10:20 -0400 Subject: [PATCH 2/3] Use context manager to open and close nc file --- city_metrix/layers/era_5_hottest_day.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/city_metrix/layers/era_5_hottest_day.py b/city_metrix/layers/era_5_hottest_day.py index 92c2982..3ec9f1c 100644 --- a/city_metrix/layers/era_5_hottest_day.py +++ b/city_metrix/layers/era_5_hottest_day.py @@ -96,14 +96,13 @@ def hourly_mean_temperature(image): }, f'download_{i}.nc') - dataarray = xr.open_dataset(f'download_{i}.nc') + with xr.open_dataset(f'download_{i}.nc') as dataarray: + # Subset times for the day + times = [valid_time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for valid_time in dataarray['valid_time'].values] + indices = [i for i, value in enumerate(times) if value in utc_times] + subset_dataarray = dataarray.isel(valid_time=indices) - # Subset times for the day - times = [valid_time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for valid_time in dataarray['valid_time'].values] - indices = [i for i, value in enumerate(times) if value in utc_times] - subset_dataarray = dataarray.isel(valid_time=indices) - - dataarray_list.append(subset_dataarray) + dataarray_list.append(subset_dataarray) # Remove local file os.remove(f'download_{i}.nc') From ca22dde6d4a7d8a64f7e69098a1251d677f65b12 Mon Sep 17 00:00:00 2001 From: Chris Rowe Date: Thu, 31 Oct 2024 20:29:27 -0400 Subject: [PATCH 3/3] Load data from file before releasing from context --- city_metrix/layers/era_5_hottest_day.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/city_metrix/layers/era_5_hottest_day.py b/city_metrix/layers/era_5_hottest_day.py index 3ec9f1c..ed4ef04 100644 --- a/city_metrix/layers/era_5_hottest_day.py +++ b/city_metrix/layers/era_5_hottest_day.py @@ -100,9 +100,9 @@ def hourly_mean_temperature(image): # Subset times for the day times = [valid_time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for valid_time in dataarray['valid_time'].values] indices = [i for i, value in enumerate(times) if value in utc_times] - subset_dataarray = dataarray.isel(valid_time=indices) + subset_dataarray = dataarray.isel(valid_time=indices).load() - dataarray_list.append(subset_dataarray) + dataarray_list.append(subset_dataarray) # Remove local file os.remove(f'download_{i}.nc')