Migrate all ODC models and prediction pipeline to Microsoft Planetary Computer
This commit is contained in:
+170
-94
@@ -1,8 +1,16 @@
|
||||
TEST_MODE = True
|
||||
RESOLUTION = 1000 if TEST_MODE else 10
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# Common imports and settings
|
||||
import os, sys
|
||||
os.environ['USE_PYGEOS'] = '0'
|
||||
os.environ["GDAL_HTTP_MAX_RETRY"] = "5"
|
||||
os.environ["GDAL_HTTP_RETRY_DELAY"] = "2"
|
||||
os.environ["GDAL_HTTP_CONNECTION_TIMEOUT"] = "10"
|
||||
os.environ["GDAL_HTTP_TIMEOUT"] = "30"
|
||||
os.environ["CPL_VSIL_CURL_ALLOWED_EXTENSIONS"] = ".tif,.tiff"
|
||||
os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "YES"
|
||||
from IPython.display import Markdown
|
||||
import pandas as pd
|
||||
pd.set_option("display.max_rows", None)
|
||||
@@ -13,15 +21,13 @@ import datacube
|
||||
from datacube.utils.rio import configure_s3_access
|
||||
from datacube.utils import masking
|
||||
from datacube.utils.cog import write_cog
|
||||
# https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools
|
||||
from dea_tools.plotting import display_map, rgb
|
||||
from dea_tools.datahandling import mostcommon_crs
|
||||
# removed deafrica_tools imports to avoid ipyleaflet error
|
||||
|
||||
# EASI defaults
|
||||
easinotebooksrepo = '/home/jovyan/easi-notebooks'
|
||||
easinotebooksrepo = '/home/x79/CSIROBoeingPhase4-Vietnam'
|
||||
if easinotebooksrepo not in sys.path: sys.path.append(easinotebooksrepo)
|
||||
from easi_tools import EasiDefaults, xarray_object_size, notebook_utils, unset_cachingproxy
|
||||
from easi_tools.load_s2l2a import load_s2l2a_with_offset
|
||||
# from easi_tools.load_s2l2a import load_s2l2a_with_offset
|
||||
from dask.distributed import progress
|
||||
|
||||
# Data tools
|
||||
@@ -31,7 +37,7 @@ from datetime import datetime
|
||||
# Datacube
|
||||
from datacube.utils import masking # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py
|
||||
from odc.algo import enum_to_bool # https://github.com/opendatacube/odc-algo/blob/main/odc/algo/_masking.py
|
||||
from odc.algo import xr_reproject # https://github.com/opendatacube/odc-algo/blob/main/odc/algo/_warp.py
|
||||
# removed xr_reproject
|
||||
from datacube.utils.geometry import GeoBox, box # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/geometry/_base.py
|
||||
|
||||
# Holoviews, Datashader and Bokeh
|
||||
@@ -83,56 +89,73 @@ import joblib
|
||||
|
||||
def load_data(dc, date_range, longtitude_range, latitude_range):
|
||||
product = 's2_l2a'
|
||||
query = {
|
||||
'product': product, # Product name
|
||||
'x': longtitude_range, # "x" axis bounds
|
||||
'y': latitude_range, # "y" axis bounds
|
||||
'time': date_range, # Any parsable date strings
|
||||
}
|
||||
native_crs = notebook_utils.mostcommon_crs(dc, query)
|
||||
print(f'Most common native CRS: {native_crs}')
|
||||
measurements = ['red', 'nir', 'scl']
|
||||
|
||||
load_params = {
|
||||
'measurements': measurements, # Selected measurement or alias names
|
||||
'output_crs': native_crs, # Target EPSG code
|
||||
'resolution': (-10, 10), # Target resolution
|
||||
'group_by': 'solar_day', # Scene grouping
|
||||
'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks
|
||||
}
|
||||
data = load_s2l2a_with_offset(
|
||||
dc,
|
||||
query | load_params # Combine the two dicts that contain our search and load parameters
|
||||
bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]
|
||||
|
||||
import pystac_client
|
||||
import planetary_computer
|
||||
import odc.stac
|
||||
|
||||
catalog = pystac_client.Client.open(
|
||||
"https://planetarycomputer.microsoft.com/api/stac/v1",
|
||||
modifier=planetary_computer.sign_inplace,
|
||||
)
|
||||
search = catalog.search(
|
||||
collections=["sentinel-2-l2a"],
|
||||
bbox=bbox,
|
||||
datetime=f"{date_range[0]}/{date_range[1]}",
|
||||
)
|
||||
items = list(search.items())
|
||||
|
||||
data = odc.stac.load(
|
||||
items,
|
||||
bands=["red", "nir", "SCL"],
|
||||
bbox=bbox,
|
||||
crs="EPSG:32648",
|
||||
resolution=RESOLUTION,
|
||||
chunks={"x": 2048, "y": 2048, "time": 1},
|
||||
groupby="solar_day"
|
||||
)
|
||||
if "SCL" in data.data_vars:
|
||||
data = data.rename({"SCL": "scl"})
|
||||
return data
|
||||
|
||||
|
||||
def mask_clean(data):
|
||||
flag_name = 'scl'
|
||||
flag_desc = masking.describe_variable_flags(data[flag_name]) # Pandas dataframe
|
||||
display(flag_desc)
|
||||
display(flag_desc.loc['qa'].values[1])
|
||||
# Create a "data quality" Mask layer
|
||||
flags_def = flag_desc.loc['qa'].values[1]
|
||||
good_pixel_flags = [flags_def[str(i)] for i in [2, 4, 5, 6]] # To pass strings to enum_to_bool()
|
||||
|
||||
# enum_to_bool calculates the pixel-wise "or" of each set of pixels given by good_pixel_flags
|
||||
# 1 = good data
|
||||
# 0 = "bad" data
|
||||
good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags)
|
||||
# For Sentinel-2 L2A SCL:
|
||||
# 2: Dark Area Pixels, 4: Vegetation, 5: Not Vegetated, 6: Water
|
||||
good_pixel_mask = data['scl'].isin([2, 4, 5, 6])
|
||||
data_layer_names = [x for x in data.data_vars if x != 'scl']
|
||||
# Apply good pixel mask to blue, green, red and nir.
|
||||
# Apply good pixel mask
|
||||
result = data[data_layer_names].where(good_pixel_mask).persist()
|
||||
return result
|
||||
|
||||
|
||||
def fill_nan(ndvi, time_split):
|
||||
if len(ndvi.time) == 0:
|
||||
return ndvi
|
||||
|
||||
# If the total time duration is less than 90 days, skip seasonal splitting
|
||||
try:
|
||||
total_days = (ndvi.time[-1] - ndvi.time[0]).dt.days.item()
|
||||
if total_days < 90:
|
||||
return ndvi.bfill(dim="time").ffill(dim="time")
|
||||
except Exception:
|
||||
pass
|
||||
|
||||
rs = []
|
||||
for times in time_split:
|
||||
tmp = ndvi.sel(time=times)
|
||||
fill_ds = tmp.sel(time=times).bfill(dim='time')
|
||||
fill_ds = fill_ds.sel(time=times).ffill(dim='time')
|
||||
rs.append(fill_ds)
|
||||
try:
|
||||
tmp = ndvi.sel(time=times)
|
||||
if len(tmp.time) == 0:
|
||||
continue
|
||||
fill_ds = tmp.bfill(dim='time').ffill(dim='time')
|
||||
rs.append(fill_ds)
|
||||
except Exception:
|
||||
continue
|
||||
|
||||
if len(rs) == 0:
|
||||
return ndvi.bfill(dim="time").ffill(dim="time")
|
||||
|
||||
merged_ndvi = xr.concat([i for i in rs], dim="time")
|
||||
fill_m = merged_ndvi.bfill(dim="time")
|
||||
fill_m = fill_m.ffill(dim="time")
|
||||
@@ -144,10 +167,49 @@ def load_train_data(train_path):
|
||||
return train
|
||||
|
||||
|
||||
def load_sen1(name_vh, name_vv):
|
||||
dsvv = rioxarray.open_rasterio(name_vv)
|
||||
dsvh = rioxarray.open_rasterio(name_vh)
|
||||
return dsvh, dsvv
|
||||
def load_sen1(bbox, time_range):
|
||||
import pystac_client
|
||||
import planetary_computer
|
||||
import odc.stac
|
||||
|
||||
# Kết nối STAC Client
|
||||
catalog = pystac_client.Client.open(
|
||||
"https://planetarycomputer.microsoft.com/api/stac/v1",
|
||||
modifier=planetary_computer.sign_inplace,
|
||||
)
|
||||
|
||||
# Tìm kiếm Items
|
||||
search = catalog.search(
|
||||
collections=["sentinel-1-rtc"],
|
||||
bbox=bbox,
|
||||
datetime=time_range,
|
||||
)
|
||||
items = list(search.items())
|
||||
|
||||
# Tải dữ liệu thành xarray Dataset
|
||||
ds_s1 = odc.stac.load(
|
||||
items,
|
||||
bands=["vv", "vh"],
|
||||
bbox=bbox,
|
||||
crs="EPSG:32648",
|
||||
resolution=RESOLUTION,
|
||||
chunks={"x": 2048, "y": 2048, "time": 1}
|
||||
)
|
||||
|
||||
# Tính giá trị trung vị theo thời gian
|
||||
ds_median = ds_s1.median(dim="time").compute()
|
||||
vv = ds_median["vv"]
|
||||
vh = ds_median["vh"]
|
||||
|
||||
# Thêm chiều 'band' để giống hệt rioxarray
|
||||
vv = vv.expand_dims(dim="band")
|
||||
vh = vh.expand_dims(dim="band")
|
||||
|
||||
# Phục hồi metadata về toạ độ
|
||||
vv = vv.rio.write_crs("EPSG:32648")
|
||||
vh = vh.rio.write_crs("EPSG:32648")
|
||||
|
||||
return vh, vv
|
||||
|
||||
|
||||
def get_data_sen1_and_sen2(train, average_ndvi, dsvh, dsvv):
|
||||
@@ -267,6 +329,10 @@ def save_model(name_file, model, metadata=None, label_encoder=None):
|
||||
|
||||
|
||||
def predict(model, data_crs, ndvi, vh, vv):
|
||||
# Unpack model if it is wrapped in a dictionary (from ModelManager)
|
||||
if isinstance(model, dict) and 'model' in model:
|
||||
model = model['model']
|
||||
|
||||
data_predict = []
|
||||
for i in range(ndvi.shape[1]):
|
||||
ndvi_tmp = ndvi.isel(y=i).values
|
||||
@@ -361,21 +427,35 @@ def save_result(result, HT_MAP):
|
||||
|
||||
def load_data_sen1(dc, date_range, coordinates):
|
||||
longtitude_range, latitude_range = coordinates
|
||||
data_sen1 = dc.load(
|
||||
product="sentinel1_grd_gamma0_10m",
|
||||
x=longtitude_range,
|
||||
y=latitude_range,
|
||||
time=date_range,
|
||||
measurements=["vv", "vh"],
|
||||
output_crs="EPSG:32648",
|
||||
resolution=(-10,10),
|
||||
dask_chunks={"x":2048, "y":2048},
|
||||
skip_broken_datasets=True,
|
||||
group_by='solar_day'
|
||||
bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]
|
||||
|
||||
import pystac_client
|
||||
import planetary_computer
|
||||
import odc.stac
|
||||
|
||||
catalog = pystac_client.Client.open(
|
||||
"https://planetarycomputer.microsoft.com/api/stac/v1",
|
||||
modifier=planetary_computer.sign_inplace,
|
||||
)
|
||||
search = catalog.search(
|
||||
collections=["sentinel-1-rtc"],
|
||||
bbox=bbox,
|
||||
datetime=f"{date_range[0]}/{date_range[1]}",
|
||||
)
|
||||
items = list(search.items())
|
||||
|
||||
data_sen1 = odc.stac.load(
|
||||
items,
|
||||
bands=["vv", "vh"],
|
||||
bbox=bbox,
|
||||
crs="EPSG:32648",
|
||||
resolution=RESOLUTION,
|
||||
chunks={"x": 2048, "y": 2048, "time": 1},
|
||||
groupby="solar_day"
|
||||
)
|
||||
|
||||
notebook_utils.heading(notebook_utils.xarray_object_size(data_sen1))
|
||||
display(data_sen1)
|
||||
# notebook_utils.heading(notebook_utils.xarray_object_size(data_sen1))
|
||||
# display(data_sen1)
|
||||
dsvh = data_sen1.vh
|
||||
dsvv = data_sen1.vv
|
||||
|
||||
@@ -387,46 +467,42 @@ def calculate_average(data, time_pattern='1M'):
|
||||
|
||||
def load_data_sen2(dc, date_range, coordinates):
|
||||
longtitude_range, latitude_range = coordinates
|
||||
product = 's2_l2a'
|
||||
query = {
|
||||
'product': product, # Product name
|
||||
'x': longtitude_range, # "x" axis bounds
|
||||
'y': latitude_range, # "y" axis bounds
|
||||
'time': date_range, # Any parsable date strings
|
||||
}
|
||||
native_crs = notebook_utils.mostcommon_crs(dc, query)
|
||||
print(f'Most common native CRS: {native_crs}')
|
||||
bbox = [longtitude_range[0], latitude_range[0], longtitude_range[1], latitude_range[1]]
|
||||
|
||||
# measurements = ['red','green', 'blue', 'nir', 'scl']
|
||||
measurements = ['red', 'nir', 'scl']
|
||||
|
||||
load_params = {
|
||||
'measurements': measurements, # Selected measurement or alias names
|
||||
'output_crs': native_crs, # Target EPSG code
|
||||
'resolution': (-10, 10), # Target resolution
|
||||
'group_by': 'solar_day', # Scene grouping
|
||||
'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks
|
||||
}
|
||||
data = load_s2l2a_with_offset(
|
||||
dc,
|
||||
query | load_params # Combine the two dicts that contain our search and load parameters
|
||||
import pystac_client
|
||||
import planetary_computer
|
||||
import odc.stac
|
||||
|
||||
catalog = pystac_client.Client.open(
|
||||
"https://planetarycomputer.microsoft.com/api/stac/v1",
|
||||
modifier=planetary_computer.sign_inplace,
|
||||
)
|
||||
search = catalog.search(
|
||||
collections=["sentinel-2-l2a"],
|
||||
bbox=bbox,
|
||||
datetime=f"{date_range[0]}/{date_range[1]}",
|
||||
)
|
||||
items = list(search.items())
|
||||
|
||||
data = odc.stac.load(
|
||||
items,
|
||||
bands=["red", "nir", "SCL"],
|
||||
bbox=bbox,
|
||||
crs="EPSG:32648",
|
||||
resolution=RESOLUTION,
|
||||
chunks={"x": 2048, "y": 2048, "time": 1},
|
||||
groupby="solar_day"
|
||||
)
|
||||
if "SCL" in data.data_vars:
|
||||
data = data.rename({"SCL": "scl"})
|
||||
return data
|
||||
|
||||
def mask_cloud(data):
|
||||
flag_name = 'scl'
|
||||
flag_desc = masking.describe_variable_flags(data[flag_name]) # Pandas dataframe
|
||||
display(flag_desc.loc['qa'].values[1])
|
||||
# Create a "data quality" Mask layer
|
||||
flags_def = flag_desc.loc['qa'].values[1]
|
||||
good_pixel_flags = [flags_def[str(i)] for i in [2, 4, 5, 6]] # To pass strings to enum_to_bool()
|
||||
|
||||
# enum_to_bool calculates the pixel-wise "or" of each set of pixels given by good_pixel_flags
|
||||
# 1 = good data
|
||||
# 0 = "bad" data
|
||||
good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags)
|
||||
# For Sentinel-2 L2A SCL:
|
||||
# 2: Dark Area Pixels, 4: Vegetation, 5: Not Vegetated, 6: Water
|
||||
good_pixel_mask = data['scl'].isin([2, 4, 5, 6])
|
||||
data_layer_names = [x for x in data.data_vars if x != 'scl']
|
||||
# Apply good pixel mask to blue, green, red and nir.
|
||||
# Apply good pixel mask
|
||||
result = data[data_layer_names].where(good_pixel_mask).persist()
|
||||
return result
|
||||
|
||||
|
||||
Reference in New Issue
Block a user