diff --git a/datacube.py b/datacube.py index 1616e690..dc22d674 100644 --- a/datacube.py +++ b/datacube.py @@ -21,20 +21,14 @@ import random from datetime import datetime, timedelta + import geopandas as gpd -import pystac -import pystac_client +import numpy as np import planetary_computer as pc +import pystac_client import stackstac import xarray as xr -import numpy as np -from shapely.geometry import box, shape, Point -import leafmap -from rasterio.errors import RasterioIOError -from rasterio.warp import transform_bounds -from getpass import getpass -import os -import json +from shapely.geometry import box STAC_API = "https://planetarycomputer.microsoft.com/api/stac/v1" S2_BANDS = ["B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12", "SCL"] @@ -74,8 +68,8 @@ def get_week(year, month, day): date = datetime(year, month, day) start_of_week = date - timedelta(days=date.weekday()) end_of_week = start_of_week + timedelta(days=6) - start_date_str = start_of_week.strftime('%Y-%m-%d') - end_date_str = end_of_week.strftime('%Y-%m-%d') + start_date_str = start_of_week.strftime("%Y-%m-%d") + end_date_str = end_of_week.strftime("%Y-%m-%d") return f"{start_date_str}/{end_date_str}" @@ -132,8 +126,17 @@ def search_sentinel2(week, aoi, cloud_cover_percentage, nodata_pixel_percentage) }, {"op": "anyinteracts", "args": [{"property": "datetime"}, week]}, {"op": "=", "args": [{"property": "collection"}, "sentinel-2-l2a"]}, - {"op": "<=", "args": [{"property": "eo:cloud_cover"}, cloud_cover_percentage]}, - {"op": "<=", "args": [{"property": "s2:nodata_pixel_percentage"}, nodata_pixel_percentage]}, + { + "op": "<=", + "args": [{"property": "eo:cloud_cover"}, cloud_cover_percentage], + }, + { + "op": "<=", + "args": [ + {"property": "s2:nodata_pixel_percentage"}, + nodata_pixel_percentage, + ], + }, ], }, ) @@ -143,24 +146,28 @@ def search_sentinel2(week, aoi, cloud_cover_percentage, nodata_pixel_percentage) s2_items_gdf = gpd.GeoDataFrame.from_features(s2_items.to_dict()) - best_nodata = (s2_items_gdf[["s2:nodata_pixel_percentage"]] - .groupby(["s2:nodata_pixel_percentage"]) - .sum() - .sort_values(by="s2:nodata_pixel_percentage", ascending=True) - .index[0]) - - best_clouds = (s2_items_gdf[["eo:cloud_cover"]] - .groupby(["eo:cloud_cover"]) - .sum() - .sort_values(by="eo:cloud_cover", ascending=True) - .index[0]) - + best_nodata = ( + s2_items_gdf[["s2:nodata_pixel_percentage"]] + .groupby(["s2:nodata_pixel_percentage"]) + .sum() + .sort_values(by="s2:nodata_pixel_percentage", ascending=True) + .index[0] + ) + + best_clouds = ( + s2_items_gdf[["eo:cloud_cover"]] + .groupby(["eo:cloud_cover"]) + .sum() + .sort_values(by="eo:cloud_cover", ascending=True) + .index[0] + ) + s2_items_gdf = s2_items_gdf[s2_items_gdf["eo:cloud_cover"] == best_clouds] # Get the item ID for the filtered Sentinel 2 dataframe containing the best cloud free scene - s2_items_gdf_datatake_id = s2_items_gdf['s2:datatake_id'] + s2_items_gdf_datatake_id = s2_items_gdf["s2:datatake_id"] for item in s2_items: - if item.properties['s2:datatake_id'] == s2_items_gdf_datatake_id[0]: + if item.properties["s2:datatake_id"] == s2_items_gdf_datatake_id[0]: s2_item = item else: continue @@ -212,11 +219,13 @@ def search_sentinel1(BBOX, catalog, week): s1_gdf = gpd.GeoDataFrame.from_features(s1_items.to_dict()) s1_gdf["overlap"] = s1_gdf.intersection(box(*BBOX)).area - state = (s1_gdf[["sat:orbit_state", "overlap"]] - .groupby(["sat:orbit_state"]) - .sum() - .sort_values(by="overlap", ascending=False) - .index[0]) + state = ( + s1_gdf[["sat:orbit_state", "overlap"]] + .groupby(["sat:orbit_state"]) + .sum() + .sort_values(by="overlap", ascending=False) + .index[0] + ) s1_gdf = s1_gdf[s1_gdf["sat:orbit_state"] == state] print("Filtered Sentinel-1 orbit state: ", s1_gdf["sat:orbit_state"].unique()) print("Number of scenes filtered by orbit state: ", len(s1_gdf)) @@ -238,9 +247,7 @@ def search_dem(BBOX, catalog, epsg): Returns: - pystac.Collection: A collection of Digital Elevation Model (DEM) items filtered by specified conditions. """ - search = catalog.search( - collections=["cop-dem-glo-30"], bbox=BBOX - ) + search = catalog.search(collections=["cop-dem-glo-30"], bbox=BBOX) dem_items = search.item_collection() print(f"Found {len(dem_items)} items") @@ -332,13 +339,17 @@ def merge_datarrays(da_sen2, da_sen1, da_dem): Returns: - xr.DataArray: Merged xarray DataArray. """ - da_merge = xr.merge([da_sen2, da_sen1, da_dem], compat='override') + da_merge = xr.merge([da_sen2, da_sen1, da_dem], compat="override") print("Merged datarray: ", da_merge) - print("Time variables (S2, merged): ", da_sen2.time.values, da_merge.time.values) # da_sen1.time.values, da_dem.time.values + print( + "Time variables (S2, merged): ", da_sen2.time.values, da_merge.time.values + ) # da_sen1.time.values, da_dem.time.values return da_merge -def process(year1, year2, aoi, resolution, cloud_cover_percentage, nodata_pixel_percentage): +def process( + year1, year2, aoi, resolution, cloud_cover_percentage, nodata_pixel_percentage +): """ Process Sentinel-2, Sentinel-1, and DEM data for a specified time range, area of interest (AOI), resolution, EPSG code, cloud cover percentage, and nodata pixel percentage. @@ -359,21 +370,28 @@ def process(year1, year2, aoi, resolution, cloud_cover_percentage, nodata_pixel_ date, YEAR, MONTH, DAY, CLOUD = get_conditions(year1, year2) week = get_week(YEAR, MONTH, DAY) - catalog, s2_items, BBOX, epsg = search_sentinel2(week, aoi, cloud_cover_percentage, nodata_pixel_percentage) + catalog, s2_items, BBOX, epsg = search_sentinel2( + week, aoi, cloud_cover_percentage, nodata_pixel_percentage + ) s1_items = search_sentinel1(BBOX, catalog, week) dem_items = search_dem(BBOX, catalog, epsg) - da_sen2, da_sen1, da_dem = make_dataarrays(s2_items, s1_items, dem_items, BBOX, resolution, epsg) + da_sen2, da_sen1, da_dem = make_dataarrays( + s2_items, s1_items, dem_items, BBOX, resolution, epsg + ) da_merge = merge_datarrays(da_sen2, da_sen1, da_dem) return da_merge + # EXAMPLE -california_tile = gpd.read_file("ca.geojson") +california_tile = gpd.read_file("ca.geojson") sample = california_tile.sample(1) aoi = sample.iloc[0].geometry cloud_cover_percentage = 50 nodata_pixel_percentage = 20 -process(2017, 2023, aoi, 10, cloud_cover_percentage, nodata_pixel_percentage) # Spatial resolution of 10 metres \ No newline at end of file +process( + 2017, 2023, aoi, 10, cloud_cover_percentage, nodata_pixel_percentage +) # Spatial resolution of 10 metres