diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index 341786b..dc1780a 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -1,8 +1,6 @@ import ee -import xarray -from dask.diagnostics import ProgressBar -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection class Albedo(Layer): """ @@ -34,19 +32,33 @@ def get_data(self, bbox): def mask_and_count_clouds(s2wc, geom): s2wc = ee.Image(s2wc) geom = ee.Geometry(geom.geometry()) - is_cloud = ee.Image(s2wc.get('cloud_mask')).gt(MAX_CLOUD_PROB).rename('is_cloud') + is_cloud = (ee.Image(s2wc.get('cloud_mask')) + .gt(MAX_CLOUD_PROB) + .rename('is_cloud') + ) + nb_cloudy_pixels = is_cloud.reduceRegion( reducer=ee.Reducer.sum().unweighted(), geometry=geom, scale=self.spatial_resolution, maxPixels=1e9 ) - return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels', - nb_cloudy_pixels.getNumber('is_cloud')).divide(10000) + mask = (s2wc + .updateMask(is_cloud.eq(0)) + .set('nb_cloudy_pixels',nb_cloudy_pixels.getNumber('is_cloud')) + .divide(10000) + ) + + return mask def mask_clouds_and_rescale(im): - clouds = ee.Image(im.get('cloud_mask')).select('probability') - return im.updateMask(clouds.lt(MAX_CLOUD_PROB)).divide(10000) + clouds = ee.Image(im.get('cloud_mask') + ).select('probability') + mask = im.updateMask(clouds + .lt(MAX_CLOUD_PROB) + ).divide(10000) + + return mask def get_masked_s2_collection(roi, start, end): criteria = (ee.Filter.And( @@ -55,23 +67,29 @@ def get_masked_s2_collection(roi, start, end): )) s2 = S2.filter(criteria) # .select('B2','B3','B4','B8','B11','B12') s2c = S2C.filter(criteria) - s2_with_clouds = (ee.Join.saveFirst('cloud_mask').apply(**{ - 'primary': ee.ImageCollection(s2), - 'secondary': ee.ImageCollection(s2c), - 'condition': ee.Filter.equals(**{'leftField': 'system:index', 'rightField': 'system:index'}) - })) + s2_with_clouds = ( + ee.Join.saveFirst('cloud_mask') + .apply(**{ + 'primary': ee.ImageCollection(s2), + 'secondary': ee.ImageCollection(s2c), + 'condition': ee.Filter.equals(**{'leftField': 'system:index', 'rightField': 'system:index'}) + }) + ) def _mcc(im): return mask_and_count_clouds(im, roi) # s2_with_clouds=ee.ImageCollection(s2_with_clouds).map(_mcc) # s2_with_clouds=s2_with_clouds.limit(image_limit,'nb_cloudy_pixels') - s2_with_clouds = ee.ImageCollection(s2_with_clouds).map( - mask_clouds_and_rescale) # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE') - return ee.ImageCollection(s2_with_clouds) + s2_with_clouds = (ee.ImageCollection(s2_with_clouds) + .map(mask_clouds_and_rescale) + ) # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE') - # calculate albedo for images + s2_with_clouds_ic = ee.ImageCollection(s2_with_clouds) + + return s2_with_clouds_ic + # calculate albedo for images # weights derived from # S. Bonafoni and A. Sekertekin, "Albedo Retrieval From Sentinel-2 by New Narrow-to-Broadband Conversion Coefficients," in IEEE Geoscience and Remote Sensing Letters, vol. 17, no. 9, pp. 1618-1622, Sept. 2020, doi: 10.1109/LGRS.2020.2967085. def calc_s2_albedo(image): @@ -89,20 +107,25 @@ def calc_s2_albedo(image): 'SWIR1': image.select('B11'), 'SWIR2': image.select('B12') } - return image.expression(S2_ALBEDO_EQN, config).double().rename('albedo') + + albedo = image.expression(S2_ALBEDO_EQN, config).double().rename('albedo') + + return albedo ## S2 MOSAIC AND ALBEDO dataset = get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date) s2_albedo = dataset.map(calc_s2_albedo) albedo_mean = s2_albedo.reduce(ee.Reducer.mean()) - data = (get_image_collection( - ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo") - .albedo_mean) + albedo_mean_ic = ee.ImageCollection(albedo_mean) + data = get_image_collection( + albedo_mean_ic, + bbox, + self.spatial_resolution, + "albedo" + ).albedo_mean if self.threshold is not None: return data.where(data < self.threshold) return data - - diff --git a/city_metrix/layers/alos_dsm.py b/city_metrix/layers/alos_dsm.py index 70000eb..c22df82 100644 --- a/city_metrix/layers/alos_dsm.py +++ b/city_metrix/layers/alos_dsm.py @@ -1,6 +1,4 @@ import ee -import xee -import xarray as xr from .layer import Layer, get_image_collection @@ -16,12 +14,19 @@ def __init__(self, spatial_resolution=30, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): - dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2") - alos_dsm = ee.ImageCollection(dataset - .filterBounds(ee.Geometry.BBox(*bbox)) - .select('DSM') - .mean() - ) - data = get_image_collection(alos_dsm, bbox, self.spatial_resolution, "ALOS DSM").DSM + alos_dsm = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2") + + alos_dsm_ic = ee.ImageCollection(alos_dsm + .filterBounds(ee.Geometry.BBox(*bbox)) + .select('DSM') + .mean() + ) + + data = get_image_collection( + alos_dsm_ic, + bbox, + self.spatial_resolution, + "ALOS DSM" + ).DSM return data diff --git a/city_metrix/layers/average_net_building_height.py b/city_metrix/layers/average_net_building_height.py index d0f49f2..11799cc 100644 --- a/city_metrix/layers/average_net_building_height.py +++ b/city_metrix/layers/average_net_building_height.py @@ -1,9 +1,7 @@ -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection + class AverageNetBuildingHeight(Layer): """ @@ -23,8 +21,13 @@ def get_data(self, bbox): # GLOBE - ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A") anbh = ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A") - data = (get_image_collection( - ee.ImageCollection(anbh), bbox, self.spatial_resolution, "average net building height") - .b1) - + + anbh_ic = ee.ImageCollection(anbh) + data = get_image_collection( + anbh_ic, + bbox, + self.spatial_resolution, + "average net building height" + ).b1 + return data diff --git a/city_metrix/layers/built_up_height.py b/city_metrix/layers/built_up_height.py index ab080f5..aef268e 100644 --- a/city_metrix/layers/built_up_height.py +++ b/city_metrix/layers/built_up_height.py @@ -1,9 +1,6 @@ -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection class BuiltUpHeight(Layer): @@ -22,8 +19,15 @@ def get_data(self, bbox): # ANBH is the average height of the built surfaces, USE THIS # AGBH is the amount of built cubic meters per surface unit in the cell # ee.ImageCollection("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_R2023A") - + built_height = ee.Image("JRC/GHSL/P2023A/GHS_BUILT_H/2018") - data = get_image_collection(ee.ImageCollection(built_height), bbox, self.spatial_resolution, "built up height") - return data.built_height + built_height_ic = ee.ImageCollection(built_height) + data = get_image_collection( + built_height_ic, + bbox, + self.spatial_resolution, + "built up height" + ).built_height + + return data diff --git a/city_metrix/layers/esa_world_cover.py b/city_metrix/layers/esa_world_cover.py index a4b0f65..c214781 100644 --- a/city_metrix/layers/esa_world_cover.py +++ b/city_metrix/layers/esa_world_cover.py @@ -1,9 +1,7 @@ -from dask.diagnostics import ProgressBar -from enum import Enum -import xarray as xr import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from enum import Enum +from .layer import Layer, get_image_collection class EsaWorldCoverClass(Enum): @@ -39,19 +37,18 @@ def __init__(self, land_cover_class=None, year=2020, spatial_resolution=10, **kw def get_data(self, bbox): if self.year == 2020: - data = get_image_collection( - ee.ImageCollection("ESA/WorldCover/v100"), - bbox, - self.spatial_resolution, - "ESA world cover" - ).Map + esa_data_ic = ee.ImageCollection("ESA/WorldCover/v100") elif self.year == 2021: - data = get_image_collection( - ee.ImageCollection("ESA/WorldCover/v200"), - bbox, - self.spatial_resolution, - "ESA world cover" - ).Map + esa_data_ic = ee.ImageCollection("ESA/WorldCover/v200") + else: + raise ValueError(f'Specified year ({self.year}) is not currently supported') + + data = get_image_collection( + esa_data_ic, + bbox, + self.spatial_resolution, + "ESA world cover" + ).Map if self.land_cover_class: data = data.where(data == self.land_cover_class.value) diff --git a/city_metrix/layers/glad_lulc.py b/city_metrix/layers/glad_lulc.py index 16d3234..11afe31 100644 --- a/city_metrix/layers/glad_lulc.py +++ b/city_metrix/layers/glad_lulc.py @@ -1,7 +1,7 @@ import xarray as xr import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection class LandCoverGlad(Layer): @@ -17,8 +17,9 @@ def __init__(self, year=2020, spatial_resolution=30, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): + lcluc_ic = ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}')) data = get_image_collection( - ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}')), + lcluc_ic, bbox, self.spatial_resolution, "GLAD Land Cover" @@ -80,7 +81,8 @@ def __init__(self, year=2020, spatial_resolution=30, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): - simplified_glad = LandCoverSimplifiedGlad(year=self.year, spatial_resolution=self.spatial_resolution).get_data(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) @@ -108,8 +110,10 @@ def __init__(self, start_year=2000, end_year=2020, spatial_resolution=30, **kwar 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) + 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)) diff --git a/city_metrix/layers/high_land_surface_temperature.py b/city_metrix/layers/high_land_surface_temperature.py index 5faa2ae..9cbf5ea 100644 --- a/city_metrix/layers/high_land_surface_temperature.py +++ b/city_metrix/layers/high_land_surface_temperature.py @@ -1,11 +1,9 @@ -from .landsat_collection_2 import LandsatCollection2 -from .land_surface_temperature import LandSurfaceTemperature -from .layer import Layer - -from shapely.geometry import box import datetime import ee +from shapely.geometry import box +from .land_surface_temperature import LandSurfaceTemperature +from .layer import Layer class HighLandSurfaceTemperature(Layer): """ @@ -28,6 +26,7 @@ def get_data(self, bbox): end_date = (hottest_date + datetime.timedelta(days=45)).strftime("%Y-%m-%d") lst = LandSurfaceTemperature(start_date, end_date, self.spatial_resolution).get_data(bbox) + lst_mean = lst.mean(dim=['x', 'y']) high_lst = lst.where(lst >= (lst_mean + self.THRESHOLD_ADD)) return high_lst @@ -36,12 +35,17 @@ def get_hottest_date(self, bbox): centroid = box(*bbox).centroid dataset = ee.ImageCollection("ECMWF/ERA5/DAILY") - AirTemperature = (dataset - .filter(ee.Filter.And( - ee.Filter.date(self.start_date, self.end_date), - ee.Filter.bounds(ee.Geometry.BBox(*bbox)))) - .select(['maximum_2m_air_temperature'], ['tasmax']) - ) + + AirTemperature = ( + dataset + .filter( + ee.Filter + .And(ee.Filter.date(self.start_date, self.end_date), + ee.Filter.bounds(ee.Geometry.BBox(*bbox)) + ) + ) + .select(['maximum_2m_air_temperature'], ['tasmax']) + ) # add date as a band to image collection def addDate(image): @@ -56,8 +60,17 @@ def addDate(image): # reduce composite to get the hottest date for centroid of ROI resolution = dataset.first().projection().nominalScale() - hottest_date = str( - ee.Number(hottest.reduceRegion(ee.Reducer.firstNonNull(), ee.Geometry.Point([centroid.x, centroid.y]), resolution).get('date')).getInfo()) + hottest_date = ( + ee.Number( + hottest.reduceRegion(ee.Reducer.firstNonNull(), + ee.Geometry.Point([centroid.x, centroid.y]), + resolution + ).get('date') + ) + .getInfo() + ) # convert to date object - return datetime.datetime.strptime(hottest_date, "%Y%m%d").date() + formated_hottest_data = datetime.datetime.strptime(str(hottest_date), "%Y%m%d").date() + + return formated_hottest_data diff --git a/city_metrix/layers/impervious_surface.py b/city_metrix/layers/impervious_surface.py index ea6772e..9cf017d 100644 --- a/city_metrix/layers/impervious_surface.py +++ b/city_metrix/layers/impervious_surface.py @@ -1,9 +1,6 @@ -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection class ImperviousSurface(Layer): @@ -18,12 +15,20 @@ def __init__(self, spatial_resolution=100, **kwargs): def get_data(self, bbox): # load impervious_surface - dataset = ee.ImageCollection(ee.Image("Tsinghua/FROM-GLC/GAIA/v10").gt(0)) # change_year_index is zero if permeable as of 2018 - imperv_surf = ee.ImageCollection(dataset - .filterBounds(ee.Geometry.BBox(*bbox)) - .select('change_year_index') - .sum() - ) - - data = get_image_collection(imperv_surf, bbox, self.spatial_resolution, "imperv surf") - return data.change_year_index + # change_year_index is zero if permeable as of 2018 + impervious_surface = ee.ImageCollection(ee.Image("Tsinghua/FROM-GLC/GAIA/v10").gt(0)) + + imperv_surf_ic = ee.ImageCollection(impervious_surface + .filterBounds(ee.Geometry.BBox(*bbox)) + .select('change_year_index') + .sum() + ) + + data = get_image_collection( + imperv_surf_ic, + bbox, + self.spatial_resolution, + "imperv surf" + ).change_year_index + + return data diff --git a/city_metrix/layers/land_surface_temperature.py b/city_metrix/layers/land_surface_temperature.py index 931cb2e..93888ce 100644 --- a/city_metrix/layers/land_surface_temperature.py +++ b/city_metrix/layers/land_surface_temperature.py @@ -1,9 +1,6 @@ -from .landsat_collection_2 import LandsatCollection2 -from .layer import Layer, get_utm_zone_epsg, get_image_collection - -from dask.diagnostics import ProgressBar import ee -import xarray + +from .layer import Layer, get_image_collection class LandSurfaceTemperature(Layer): """ @@ -22,6 +19,7 @@ def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resol def get_data(self, bbox): def cloud_mask(image): qa = image.select('QA_PIXEL') + mask = qa.bitwiseAnd(1 << 3).Or(qa.bitwiseAnd(1 << 4)) return image.updateMask(mask.Not()) @@ -30,13 +28,22 @@ def apply_scale_factors(image): return thermal_band l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") - l8_st = l8 \ - .select('ST_B10', 'QA_PIXEL') \ - .filter(ee.Filter.date(self.start_date, self.end_date)) \ - .filterBounds(ee.Geometry.BBox(*bbox)) \ - .map(cloud_mask) \ - .map(apply_scale_factors) \ - .reduce(ee.Reducer.mean()) - - data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.spatial_resolution, "LST").ST_B10_mean + + l8_st = (l8 + .select('ST_B10', 'QA_PIXEL') + .filter(ee.Filter.date(self.start_date, self.end_date)) + .filterBounds(ee.Geometry.BBox(*bbox)) + .map(cloud_mask) + .map(apply_scale_factors) + .reduce(ee.Reducer.mean()) + ) + + l8_st_ic = ee.ImageCollection(l8_st) + data = get_image_collection( + l8_st_ic, + bbox, + self.spatial_resolution, + "LST" + ).ST_B10_mean + return data diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index cad46bf..9beae31 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -2,7 +2,6 @@ from abc import abstractmethod from typing import Union, Tuple from uuid import uuid4 -from osgeo import gdal import ee import boto3 diff --git a/city_metrix/layers/nasa_dem.py b/city_metrix/layers/nasa_dem.py index b5ac45d..d3840d3 100644 --- a/city_metrix/layers/nasa_dem.py +++ b/city_metrix/layers/nasa_dem.py @@ -1,6 +1,4 @@ import ee -import xee -import xarray as xr from .layer import Layer, get_image_collection @@ -16,12 +14,20 @@ def __init__(self, spatial_resolution=30, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): - dataset = ee.Image("NASA/NASADEM_HGT/001") - nasa_dem = ee.ImageCollection(ee.ImageCollection(dataset) - .filterBounds(ee.Geometry.BBox(*bbox)) - .select('elevation') - .mean() - ) - data = get_image_collection(nasa_dem, bbox, self.spatial_resolution, "NASA DEM").elevation + nasa_dem = ee.Image("NASA/NASADEM_HGT/001") + + nasa_dem_elev = (ee.ImageCollection(nasa_dem) + .filterBounds(ee.Geometry.BBox(*bbox)) + .select('elevation') + .mean() + ) + + nasa_dem_elev_ic = ee.ImageCollection(nasa_dem_elev) + data = get_image_collection( + nasa_dem_elev_ic, + bbox, + self.spatial_resolution, + "NASA DEM" + ).elevation return data diff --git a/city_metrix/layers/natural_areas.py b/city_metrix/layers/natural_areas.py index 5efe4a8..fdf499b 100644 --- a/city_metrix/layers/natural_areas.py +++ b/city_metrix/layers/natural_areas.py @@ -1,4 +1,3 @@ -import xarray as xr from xrspatial.classify import reclassify from .layer import Layer @@ -32,7 +31,11 @@ def get_data(self, bbox): } # Perform the reclassification - reclassified_data = reclassify(esa_world_cover, bins=list(reclass_map.keys()), new_values=list(reclass_map.values())) + reclassified_data = reclassify( + esa_world_cover, + bins=list(reclass_map.keys()), + new_values=list(reclass_map.values()) + ) # Apply the original CRS (Coordinate Reference System) reclassified_data = reclassified_data.rio.write_crs(esa_world_cover.rio.crs, inplace=True) diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py index fce3021..b31021d 100644 --- a/city_metrix/layers/ndvi_sentinel2_gee.py +++ b/city_metrix/layers/ndvi_sentinel2_gee.py @@ -1,4 +1,5 @@ import ee + from .layer import Layer, get_image_collection class NdviSentinel2(Layer): @@ -32,6 +33,7 @@ def calculate_ndvi(image): return image.addBands(ndvi) s2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED") + ndvi = (s2 .filterBounds(ee.Geometry.BBox(*bbox)) .filterDate(start_date, end_date) @@ -41,8 +43,11 @@ def calculate_ndvi(image): ndvi_mosaic = ndvi.qualityMosaic('NDVI') - ic = ee.ImageCollection(ndvi_mosaic) - ndvi_data = (get_image_collection(ic, bbox, self.spatial_resolution, "NDVI") - .NDVI) + ndvi_mosaic_ic = ee.ImageCollection(ndvi_mosaic) + ndvi_data = get_image_collection( + ndvi_mosaic_ic, + bbox, + self.spatial_resolution, "NDVI" + ).NDVI return ndvi_data diff --git a/city_metrix/layers/smart_surface_lulc.py b/city_metrix/layers/smart_surface_lulc.py index e4759a6..2e0f083 100644 --- a/city_metrix/layers/smart_surface_lulc.py +++ b/city_metrix/layers/smart_surface_lulc.py @@ -1,9 +1,7 @@ import xarray as xr import numpy as np import pandas as pd -import geopandas as gpd from shapely.geometry import CAP_STYLE, JOIN_STYLE -from shapely.geometry import box from exactextract import exact_extract from geocube.api.core import make_geocube import warnings @@ -20,15 +18,16 @@ class SmartSurfaceLULC(Layer): - def __init__(self, land_cover_class=None, **kwargs): + def __init__(self, land_cover_class=None, spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.land_cover_class = land_cover_class + self.spatial_resolution = spatial_resolution def get_data(self, bbox): crs = get_utm_zone_epsg(bbox) # ESA world cover - esa_world_cover = EsaWorldCover(year=2021).get_data(bbox) + esa_world_cover = EsaWorldCover(year=2021, spatial_resolution = self.spatial_resolution).get_data(bbox) # ESA reclass and upsample def get_data_esa_reclass(esa_world_cover): reclass_map = { @@ -46,7 +45,11 @@ def get_data_esa_reclass(esa_world_cover): } # Perform the reclassification - reclassified_esa = reclassify(esa_world_cover, bins=list(reclass_map.keys()), new_values=list(reclass_map.values())) + reclassified_esa = reclassify( + esa_world_cover, + bins=list(reclass_map.keys()), + new_values=list(reclass_map.values()) + ) esa_1m = reclassified_esa.rio.reproject( dst_crs=crs, @@ -74,7 +77,8 @@ def get_data_esa_reclass(esa_world_cover): if len(roads_osm) > 0: roads_osm['lanes'] = pd.to_numeric(roads_osm['lanes'], errors='coerce') # Get the average number of lanes per highway class - lanes = (roads_osm.drop(columns='geometry') + lanes = (roads_osm + .drop(columns='geometry') .groupby('highway') # Calculate average and round up .agg(avg_lanes=('lanes', lambda x: np.ceil(np.nanmean(x)) if not np.isnan(x).all() else np.NaN)) @@ -93,11 +97,16 @@ def get_data_esa_reclass(esa_world_cover): # https://nacto.org/publication/urban-street-design-guide/street-design-elements/lane-width/#:~:text=wider%20lane%20widths.-,Lane%20widths%20of%2010%20feet%20are%20appropriate%20in%20urban%20areas,be%20used%20in%20each%20direction # cap is flat to the terminus of the road # join style is mitred so intersections are squared - roads_osm['geometry'] = roads_osm.apply(lambda row: row['geometry'].buffer( - row['lanes'] * 3.048 / 2, - cap_style=CAP_STYLE.flat, - join_style=JOIN_STYLE.mitre), - axis=1 + roads_osm['geometry'] = ( + roads_osm + .apply( + lambda row: row['geometry'] + .buffer( + row['lanes'] * 3.048 / 2, + cap_style=CAP_STYLE.flat, + join_style=JOIN_STYLE.mitre + ),axis=1 + ) ) else: # Add value field (30) @@ -179,14 +188,20 @@ def classify_building(building): if len(buildings) > 0: buildings['Value'] = buildings.apply(classify_building, axis=1) - # Parking parking_osm = OpenStreetMap(osm_class=OpenStreetMapClass.PARKING).get_data(bbox).reset_index() parking_osm['Value'] = 50 - # combine features: open space, water, road, building, parking - feature_df = pd.concat([open_space_osm[['geometry', 'Value']], water_osm[['geometry', 'Value']], roads_osm[['geometry', 'Value']], buildings[['geometry', 'Value']], parking_osm[['geometry', 'Value']]], axis=0) + feature_df = pd.concat( + [open_space_osm[['geometry', 'Value']], + water_osm[['geometry', 'Value']], + roads_osm[['geometry', 'Value']], + buildings[['geometry', 'Value']], + parking_osm[['geometry', 'Value']] + ], axis=0 + ) + # rasterize if feature_df.empty: feature_1m = xr.zeros_like(esa_1m) diff --git a/city_metrix/layers/tree_canopy_height.py b/city_metrix/layers/tree_canopy_height.py index fee1697..f499552 100644 --- a/city_metrix/layers/tree_canopy_height.py +++ b/city_metrix/layers/tree_canopy_height.py @@ -1,10 +1,7 @@ -from .layer import Layer, get_utm_zone_epsg, get_image_collection - -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee +from .layer import Layer, get_image_collection + class TreeCanopyHeight(Layer): """ Attributes: @@ -20,10 +17,19 @@ def __init__(self, spatial_resolution=1, **kwargs): def get_data(self, bbox): canopy_ht = ee.ImageCollection("projects/meta-forest-monitoring-okw37/assets/CanopyHeight") + # aggregate time series into a single image - canopy_ht = canopy_ht.reduce(ee.Reducer.mean()).rename("cover_code") + canopy_ht_img = (canopy_ht + .reduce(ee.Reducer.mean()) + .rename("cover_code") + ) - data = get_image_collection(ee.ImageCollection(canopy_ht), bbox, - self.spatial_resolution, "tree canopy height") + canopy_ht_ic = ee.ImageCollection(canopy_ht_img) + data = get_image_collection( + canopy_ht_ic, + bbox, + self.spatial_resolution, + "tree canopy height" + ).cover_code - return data.cover_code + return data diff --git a/city_metrix/layers/tree_cover.py b/city_metrix/layers/tree_cover.py index 98bc481..52fc8d6 100644 --- a/city_metrix/layers/tree_cover.py +++ b/city_metrix/layers/tree_cover.py @@ -1,10 +1,7 @@ -from .layer import Layer, get_utm_zone_epsg, get_image_collection - -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee +from .layer import Layer, get_image_collection + class TreeCover(Layer): """ Merged tropical and nontropical tree cover from WRI @@ -24,11 +21,21 @@ def __init__(self, min_tree_cover=None, max_tree_cover=None, spatial_resolution= def get_data(self, bbox): tropics = ee.ImageCollection('projects/wri-datalab/TropicalTreeCover') - nontropics = ee.ImageCollection('projects/wri-datalab/TTC-nontropics') - merged_ttc = tropics.merge(nontropics) - ttc_image = merged_ttc.reduce(ee.Reducer.mean()).rename('ttc') - - data = get_image_collection(ee.ImageCollection(ttc_image), bbox, self.spatial_resolution, "tree cover").ttc + non_tropics = ee.ImageCollection('projects/wri-datalab/TTC-nontropics') + + merged_ttc = tropics.merge(non_tropics) + ttc_image = (merged_ttc + .reduce(ee.Reducer.mean()) + .rename('ttc') + ) + + ttc_ic = ee.ImageCollection(ttc_image) + data = get_image_collection( + ttc_ic, + bbox, + self.spatial_resolution, + "tree cover" + ).ttc if self.min_tree_cover is not None: data = data.where(data >= self.min_tree_cover) @@ -36,4 +43,3 @@ def get_data(self, bbox): data = data.where(data <= self.max_tree_cover) return data - diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index fe69c75..02c737c 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -1,9 +1,6 @@ -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection class UrbanLandUse(Layer): @@ -19,21 +16,28 @@ def __init__(self, band='lulc', spatial_resolution=5, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): - dataset = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") + ulu = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") + # ImageCollection didn't cover the globe - if dataset.filterBounds(ee.Geometry.BBox(*bbox)).size().getInfo() == 0: - ulu = ee.ImageCollection(ee.Image.constant(0) + if ulu.filterBounds(ee.Geometry.BBox(*bbox)).size().getInfo() == 0: + ulu_ic = ee.ImageCollection(ee.Image + .constant(0) .clip(ee.Geometry.BBox(*bbox)) .rename('lulc') ) else: - ulu = ee.ImageCollection(dataset + ulu_ic = ee.ImageCollection(ulu .filterBounds(ee.Geometry.BBox(*bbox)) .select(self.band) .reduce(ee.Reducer.firstNonNull()) .rename('lulc') ) - data = get_image_collection(ulu, bbox, self.spatial_resolution, "urban land use").lulc + data = get_image_collection( + ulu_ic, + bbox, + self.spatial_resolution, + "urban land use" + ).lulc return data diff --git a/city_metrix/layers/world_pop.py b/city_metrix/layers/world_pop.py index 30fb6d8..700010a 100644 --- a/city_metrix/layers/world_pop.py +++ b/city_metrix/layers/world_pop.py @@ -1,9 +1,6 @@ -from dask.diagnostics import ProgressBar -import xarray as xr -import xee import ee -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_image_collection class WorldPop(Layer): @@ -27,23 +24,44 @@ def get_data(self, bbox): if not self.agesex_classes: # total population dataset = ee.ImageCollection('WorldPop/GP/100m/pop') - world_pop = ee.ImageCollection(dataset - .filterBounds(ee.Geometry.BBox(*bbox)) - .filter(ee.Filter.inList('year', [self.year])) - .select('population') - .mean() - ) - data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop").population + world_pop_ic = ee.ImageCollection( + dataset + .filterBounds(ee.Geometry.BBox(*bbox)) + .filter(ee.Filter.inList('year', [self.year])) + .select('population') + .mean() + ) + + data = get_image_collection( + world_pop_ic, + bbox, + self.spatial_resolution, + "world pop" + ).population else: # sum population for selected age-sex groups - dataset = ee.ImageCollection('WorldPop/GP/100m/pop_age_sex') - world_pop = dataset.filterBounds(ee.Geometry.BBox(*bbox))\ - .filter(ee.Filter.inList('year', [self.year]))\ - .select(self.agesex_classes)\ - .mean() - world_pop = ee.ImageCollection(world_pop.reduce(ee.Reducer.sum()).rename('sum_age_sex_group')) - data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop age sex").sum_age_sex_group + world_pop_age_sex = ee.ImageCollection('WorldPop/GP/100m/pop_age_sex') + + world_pop_age_sex_year = (world_pop_age_sex + .filterBounds(ee.Geometry.BBox(*bbox)) + .filter(ee.Filter.inList('year', [self.year])) + .select(self.agesex_classes) + .mean() + ) + + world_pop_group_ic = ee.ImageCollection( + world_pop_age_sex_year + .reduce(ee.Reducer.sum()) + .rename('sum_age_sex_group') + ) + + data = get_image_collection( + world_pop_group_ic, + bbox, + self.spatial_resolution, + "world pop age sex" + ).sum_age_sex_group return data diff --git a/environment.yml b/environment.yml index 0350113..3953796 100644 --- a/environment.yml +++ b/environment.yml @@ -6,6 +6,7 @@ dependencies: - earthengine-api=0.1.411 - geocube=0.4.2 - geopandas=0.14.4 + - xarray=2024.7.0 - rioxarray=0.15.0 - odc-stac=0.3.8 - pystac-client=0.7.5 diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py index d72876d..6b33b51 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py @@ -24,18 +24,13 @@ CUSTOM_DUMP_DIRECTORY = None def pytest_configure(config): - qgis_project_file = 'layers_for_br_lauro_de_freitas.qgz' if RUN_DUMPS is True: source_folder = os.path.dirname(__file__) target_folder = get_target_folder_path() - create_target_folder(target_folder, True) + create_target_folder(target_folder, False) - source_qgis_file = os.path.join(source_folder, qgis_project_file) - target_qgis_file = os.path.join(target_folder, qgis_project_file) - shutil.copyfile(source_qgis_file, target_qgis_file) - - print("\n\033[93m QGIS project file and layer files written to folder %s.\033[0m\n" % target_folder) + print("\n\033[93m Layer files written to folder %s.\033[0m\n" % target_folder) @pytest.fixture def target_folder(): diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz b/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz deleted file mode 100644 index 1cc53ce..0000000 Binary files a/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz and /dev/null differ diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_all_layers.py similarity index 81% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py rename to tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_all_layers.py index 9c38211..2be8c4d 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_all_layers.py @@ -2,27 +2,7 @@ # Execution configuration is specified in the conftest file import pytest -from city_metrix.layers import ( - Albedo, - AlosDSM, - AverageNetBuildingHeight, - EsaWorldCover, - HighLandSurfaceTemperature, - LandsatCollection2, - LandSurfaceTemperature, - NasaDEM, - NaturalAreas, - OpenBuildings, - OpenStreetMap, - OvertureBuildings, - Sentinel2Level2, - NdviSentinel2, - SmartSurfaceLULC, - TreeCanopyHeight, - TreeCover, - UrbanLandUse, - WorldPop, Layer, ImperviousSurface -) +from city_metrix.layers import * from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated, get_file_count_in_folder from ...tools.general_tools import get_class_default_spatial_resolution @@ -34,30 +14,6 @@ def test_write_albedo(target_folder, bbox_info, target_spatial_resolution_multip Albedo(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) -@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_albedo_tiled_unbuffered(target_folder, bbox_info, target_spatial_resolution_multiplier): - file_path = prep_output_path(target_folder, 'albedo_tiled.tif') - target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(Albedo()) - (Albedo(spatial_resolution=target_resolution). - write(bbox_info.bounds, file_path, tile_degrees=0.01, buffer_size=None)) - file_count = get_file_count_in_folder(file_path) - - expected_file_count = 5 # includes 4 tiles and one geojson file - assert file_count == expected_file_count - -@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_albedo_tiled_buffered(target_folder, bbox_info, target_spatial_resolution_multiplier): - buffer_degrees = 0.001 - file_path = prep_output_path(target_folder, 'albedo_tiled_buffered.tif') - target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(Albedo()) - (Albedo(spatial_resolution=target_resolution). - write(bbox_info.bounds, file_path, tile_degrees=0.01, buffer_size=buffer_degrees)) - file_count = get_file_count_in_folder(file_path) - - expected_file_count = 6 # includes 4 tiles and two geojson files - assert file_count == expected_file_count - - @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') def test_write_alos_dsm(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'alos_dsm.tif') @@ -72,6 +28,13 @@ def test_write_average_net_building_height(target_folder, bbox_info, target_spat AverageNetBuildingHeight(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_built_up_height(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, 'built_up_height.tif') + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(BuiltUpHeight()) + BuiltUpHeight(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') def test_write_esa_world_cover(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'esa_world_cover.tif') @@ -79,6 +42,13 @@ def test_write_esa_world_cover(target_folder, bbox_info, target_spatial_resoluti EsaWorldCover(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_glad_lulc(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, 'glad_lulc.tif') + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(LandCoverGlad()) + LandCoverGlad(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') def test_write_high_land_surface_temperature(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'high_land_surface_temperature.tif') @@ -90,14 +60,7 @@ def test_write_high_land_surface_temperature(target_folder, bbox_info, target_sp def test_write_impervious_surface(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'impervious_surface.tif') target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(ImperviousSurface()) - LandSurfaceTemperature(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) - assert verify_file_is_populated(file_path) - -@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_land_surface_temperature(target_folder, bbox_info, target_spatial_resolution_multiplier): - file_path = prep_output_path(target_folder, 'land_surface_temperature.tif') - target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(LandSurfaceTemperature()) - LandSurfaceTemperature(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) + ImperviousSurface(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) # TODO Class is no longer used, but may be useful later @@ -158,9 +121,9 @@ def test_write_overture_buildings(target_folder, bbox_info, target_spatial_resol @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') def test_write_smart_surface_lulc(target_folder, bbox_info, target_spatial_resolution_multiplier): - # Note: spatial_resolution not implemented for this raster class file_path = prep_output_path(target_folder, 'smart_surface_lulc.tif') - SmartSurfaceLULC().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(SmartSurfaceLULC()) + SmartSurfaceLULC(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_other.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_other.py new file mode 100644 index 0000000..d4b0a23 --- /dev/null +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_other.py @@ -0,0 +1,31 @@ +# This code is mostly intended for manual execution +# Execution configuration is specified in the conftest file +import pytest + +from city_metrix.layers import * +from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated, get_file_count_in_folder +from ...tools.general_tools import get_class_default_spatial_resolution + + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_albedo_tiled_unbuffered(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, 'albedo_tiled.tif') + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(Albedo()) + (Albedo(spatial_resolution=target_resolution). + write(bbox_info.bounds, file_path, tile_degrees=0.01, buffer_size=None)) + file_count = get_file_count_in_folder(file_path) + + expected_file_count = 5 # includes 4 tiles and one geojson file + assert file_count == expected_file_count + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_albedo_tiled_buffered(target_folder, bbox_info, target_spatial_resolution_multiplier): + buffer_degrees = 0.001 + file_path = prep_output_path(target_folder, 'albedo_tiled_buffered.tif') + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(Albedo()) + (Albedo(spatial_resolution=target_resolution). + write(bbox_info.bounds, file_path, tile_degrees=0.01, buffer_size=buffer_degrees)) + file_count = get_file_count_in_folder(file_path) + + expected_file_count = 6 # includes 4 tiles and two geojson files + assert file_count == expected_file_count \ No newline at end of file diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py new file mode 100644 index 0000000..8964b86 --- /dev/null +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py @@ -0,0 +1,111 @@ +# This code is mostly intended for manual execution +# Execution configuration is specified in the conftest file +import pytest + +from city_metrix.layers import * +from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated, get_file_count_in_folder + +TARGET_RESOLUTION = 5 + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_albedo_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'albedo_{TARGET_RESOLUTION}m.tif') + Albedo(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_alos_dsm_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'alos_dsm_{TARGET_RESOLUTION}m.tif') + AlosDSM(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_average_net_building_height_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'average_net_building_height_{TARGET_RESOLUTION}m.tif') + AverageNetBuildingHeight(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_built_up_height_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'built_up_height_{TARGET_RESOLUTION}m.tif') + BuiltUpHeight(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_esa_world_cover_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'esa_world_cover_{TARGET_RESOLUTION}m.tif') + EsaWorldCover(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_glad_lulc_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'glad_lulc_{TARGET_RESOLUTION}m.tif') + LandCoverGlad(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_high_land_surface_temperature_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'high_land_surface_temperature_{TARGET_RESOLUTION}m.tif') + HighLandSurfaceTemperature(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_impervious_surface_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'impervious_surface_{TARGET_RESOLUTION}m.tif') + ImperviousSurface(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_land_surface_temperature_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'land_surface_temperature_{TARGET_RESOLUTION}m.tif') + LandSurfaceTemperature(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_nasa_dem_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'nasa_dem_{TARGET_RESOLUTION}m.tif') + NasaDEM(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_natural_areas_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'natural_areas_{TARGET_RESOLUTION}m.tif') + NaturalAreas(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_ndvi_sentinel2_gee_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'ndvi_sentinel2_gee_{TARGET_RESOLUTION}m.tif') + NdviSentinel2(year=2023, spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_smart_surface_lulc_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'smart_surface_lulc_{TARGET_RESOLUTION}m.tif') + SmartSurfaceLULC(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_tree_canopy_height_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'tree_canopy_height_{TARGET_RESOLUTION}m.tif') + TreeCanopyHeight(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_tree_cover_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'tree_cover_{TARGET_RESOLUTION}m.tif') + TreeCover(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_urban_land_use_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'urban_land_use_{TARGET_RESOLUTION}m.tif') + UrbanLandUse(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_world_pop_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'world_pop_{TARGET_RESOLUTION}m.tif') + WorldPop(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None) + assert verify_file_is_populated(file_path) + diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 3fe415b..359bb28 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,4 +1,5 @@ import pytest +import xarray as xr import numpy as np from skimage.metrics import structural_similarity as ssim from pyproj import CRS @@ -286,7 +287,10 @@ def _get_sample_data(class_instance, bbox, downsize_factor): return default_res_data, downsized_res_data, downsized_resolution, estimated_actual_resolution def _get_crs_from_image_data(image_data): - crs_string = image_data.rio.crs.data['init'] + if isinstance(image_data, xr.DataArray): + crs_string = image_data.rio.crs.wkt + else: + crs_string = image_data.crs crs = CRS.from_string(crs_string) return crs