diff --git a/.github/requirements.txt b/.github/requirements.txt index 90f10fe..f5069a9 100644 --- a/.github/requirements.txt +++ b/.github/requirements.txt @@ -1,6 +1,6 @@ earthengine-api==0.1.408 geocube==0.4.2 -geopandas==0.14.4 +geopandas==0.14.1 rioxarray==0.15.0 odc-stac==0.3.8 pystac-client==0.7.5 @@ -15,7 +15,7 @@ s3fs==2024.5.0 geemap==0.32.0 pip==23.3.1 boto3==1.34.124 -scikit-learn==1.5.0 +scikit-learn==1.5.1 scikit-image==0.24.0 overturemaps==0.6.0 exactextract==0.2.0.dev252 diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py index 669e727..7e5e19e 100644 --- a/city_metrix/layers/__init__.py +++ b/city_metrix/layers/__init__.py @@ -19,3 +19,4 @@ from .alos_dsm import AlosDSM from .overture_buildings import OvertureBuildings from .nasa_dem import NasaDEM +from .impervious_surface import ImperviousSurface diff --git a/city_metrix/layers/impervious_surface.py b/city_metrix/layers/impervious_surface.py new file mode 100644 index 0000000..74e624d --- /dev/null +++ b/city_metrix/layers/impervious_surface.py @@ -0,0 +1,23 @@ +from dask.diagnostics import ProgressBar +import xarray as xr +import xee +import ee + +from .layer import Layer, get_utm_zone_epsg, get_image_collection + + +class ImperviousSurface(Layer): + def __init__(self, **kwargs): + super().__init__(**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, 100, "imperv surf") + return data.change_year_index diff --git a/city_metrix/layers/smart_surface_lulc.py b/city_metrix/layers/smart_surface_lulc.py index 2ab2a02..1f39496 100644 --- a/city_metrix/layers/smart_surface_lulc.py +++ b/city_metrix/layers/smart_surface_lulc.py @@ -13,6 +13,7 @@ from .layer import Layer, get_utm_zone_epsg, create_fishnet_grid, MAX_TILE_SIZE from .open_street_map import OpenStreetMap, OpenStreetMapClass +from .overture_buildings import OvertureBuildings from ..models.building_classifier.building_classifier import BuildingClassifier @@ -33,12 +34,12 @@ def get_data(self, bbox): # Open space open_space_osm = OpenStreetMap(osm_class=OpenStreetMapClass.OPEN_SPACE_HEAT).get_data(bbox).to_crs(crs).reset_index() - open_space_osm['Value'] = np.int8(10) + open_space_osm['Value'] = 10 # Water water_osm = OpenStreetMap(osm_class=OpenStreetMapClass.WATER).get_data(bbox).to_crs(crs).reset_index() - water_osm['Value'] = np.int8(20) + water_osm['Value'] = 20 # Roads @@ -59,7 +60,7 @@ def get_data(self, bbox): roads_osm['lanes'] = roads_osm['lanes'].fillna(roads_osm['avg_lanes']) # Add value field (30) - roads_osm['Value'] = np.int8(30) + roads_osm['Value'] = 30 # Buffer roads by lanes * 10 ft (3.048 m) # 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 @@ -73,14 +74,14 @@ def get_data(self, bbox): ) else: # Add value field (30) - roads_osm['Value'] = np.int8(30) + roads_osm['Value'] = 30 # Building ulu_lulc_1m = BuildingClassifier().get_data_ulu(bbox, crs, esa_1m) anbh_1m = BuildingClassifier().get_data_anbh(bbox, esa_1m) # get building features - buildings = BuildingClassifier().get_data_buildings(bbox, crs) + buildings = OvertureBuildings().get_data(bbox).to_crs(crs) # extract ULU, ANBH, and Area_m buildings['ULU'] = exact_extract(ulu_lulc_1m, buildings, ["majority"], output='pandas')['majority'] buildings['ANBH'] = exact_extract(anbh_1m, buildings, ["mean"], output='pandas')['mean'] @@ -112,7 +113,7 @@ def get_data(self, bbox): # Parking parking_osm = OpenStreetMap(osm_class=OpenStreetMapClass.PARKING).get_data(bbox).to_crs(crs).reset_index() - parking_osm['Value'] = np.int8(50) + parking_osm['Value'] = 50 # combine features: open space, water, road, building, parking diff --git a/city_metrix/models/building_classifier/building_classifier.pkl b/city_metrix/models/building_classifier/building_classifier.pkl index d1cee21..4907a22 100644 Binary files a/city_metrix/models/building_classifier/building_classifier.pkl and b/city_metrix/models/building_classifier/building_classifier.pkl differ diff --git a/city_metrix/models/building_classifier/building_classifier.py b/city_metrix/models/building_classifier/building_classifier.py index 2241509..816b690 100644 --- a/city_metrix/models/building_classifier/building_classifier.py +++ b/city_metrix/models/building_classifier/building_classifier.py @@ -16,8 +16,6 @@ from ...layers.esa_world_cover import EsaWorldCover, EsaWorldCoverClass from ...layers.urban_land_use import UrbanLandUse from ...layers.average_net_building_height import AverageNetBuildingHeight -from ...layers.open_street_map import OpenStreetMap, OpenStreetMapClass -from ...layers.open_buildings import OpenBuildings class BuildingClassifier(Layer): @@ -52,8 +50,8 @@ def get_data_esa_reclass(self, bbox, crs): # Perform the reclassification reclassified_esa = reclassify(esa_world_cover, bins=list(reclass_map.keys()), new_values=list(reclass_map.values())) - # Convert to int8 and chunk the data for Dask processing - reclassified_esa = reclassified_esa.astype(np.int8).chunk({'x': 512, 'y': 512}) + # Chunk the data for Dask processing + reclassified_esa = reclassified_esa.chunk({'x': 512, 'y': 512}) reclassified_esa = reclassified_esa.rio.write_crs(esa_world_cover.rio.crs, inplace=True) @@ -81,8 +79,8 @@ def get_data_ulu(self, bbox, crs, snap_to): for from_val, to_val in mapping.items(): ulu_lulc = ulu_lulc.where(ulu_lulc != from_val, to_val) - # Convert to int8 and chunk the data for Dask processing - ulu_lulc = ulu_lulc.astype(np.int8).chunk({'x': 512, 'y': 512}) + # Chunk the data for Dask processing + ulu_lulc = ulu_lulc.chunk({'x': 512, 'y': 512}) ####### 1-Non-residential as default # 0-Unclassified as nodata @@ -111,28 +109,10 @@ def get_data_anbh(self, bbox, snap_to): return anbh_1m - def get_data_buildings(self, bbox, crs): - # OSM buildings - building_osm = OpenStreetMap(osm_class=OpenStreetMapClass.BUILDING).get_data(bbox).to_crs(crs).reset_index(drop=True) - # Google-Microsoft Open Buildings Dataset buildings - openbuilds = OpenBuildings(country='USA').get_data(bbox).to_crs(crs).reset_index(drop=True) - - # Intersect buildings and keep the open buildings that don't intersect OSM buildings - intersect_buildings = gpd.sjoin(building_osm, openbuilds, how='inner', predicate='intersects') - openbuilds_non_intersect = openbuilds.loc[~openbuilds.index.isin(intersect_buildings.index)] - - buildings = pd.concat([building_osm['geometry'], openbuilds_non_intersect['geometry']], ignore_index=True).reset_index() - # Get rid of any 3d geometries that cause a problem - buildings = buildings[~buildings['geometry'].apply(lambda geom: 'Z' in geom.geom_type)] - - # Value not start with 0 - buildings['Value'] = buildings['index'] + 1 - - return buildings def rasterize_polygon(self, gdf, snap_to): if gdf.empty: - raster = np.full(snap_to.shape, 0, dtype=np.int8) + raster = np.full(snap_to.shape, 0) raster = xr.DataArray(raster, dims=snap_to.dims, coords=snap_to.coords) return raster.rio.write_crs(snap_to.rio.crs, inplace=True) @@ -141,7 +121,7 @@ def rasterize_polygon(self, gdf, snap_to): vector_data=gdf, measurements=["Value"], like=snap_to, - fill=np.int8(0) + fill=0 ).Value return raster.rio.reproject_match(snap_to) @@ -169,7 +149,7 @@ def building_classifier_tree(self): # set classifier parameters clf = DecisionTreeClassifier(max_depth=5) # encode labels - buildings_sample['Slope_encoded'] = buildings_sample['Slope'].map({'low': np.int8(42), 'high': np.int8(40)}) + buildings_sample['Slope_encoded'] = buildings_sample['Slope'].map({'low': 42, 'high': 40}) # Select these rows for the training set build_train = buildings_sample[buildings_sample['Model']=='training'] diff --git a/environment.yml b/environment.yml index 18ac340..a064834 100644 --- a/environment.yml +++ b/environment.yml @@ -21,7 +21,7 @@ dependencies: - geemap=0.32.0 - pip=23.3.1 - boto3=1.34.124 - - scikit-learn=1.5.0 + - scikit-learn=1.5.1 - scikit-image==0.24.0 - exactextract=0.2.0.dev252 - pip: diff --git a/setup.py b/setup.py index ff88cbd..7d9242c 100644 --- a/setup.py +++ b/setup.py @@ -30,7 +30,7 @@ "boto3", "exactextract<=0.2.0.dev252", "overturemaps", - "scikit-learn>=1.5.0", + "scikit-learn>=1.5.1", "scikit-image>=0.24.0" ], ) diff --git a/tests/test_layers.py b/tests/test_layers.py index bfcf4a0..27cb7c5 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -9,6 +9,7 @@ EsaWorldCover, EsaWorldCoverClass, HighLandSurfaceTemperature, + ImperviousSurface, LandsatCollection2, LandSurfaceTemperature, NasaDEM, @@ -82,6 +83,11 @@ def test_high_land_surface_temperature(): assert data.any() +def test_impervious_surface(): + data = ImperviousSurface().get_data(BBOX) + assert data.any() + + def test_land_surface_temperature(): mean_lst = LandSurfaceTemperature().get_data(BBOX).mean() assert mean_lst