diff --git a/api/app/auto_spatial_advisory/process_hfi.py b/api/app/auto_spatial_advisory/process_hfi.py index 5b8cfcbb2..1c9eceaa1 100644 --- a/api/app/auto_spatial_advisory/process_hfi.py +++ b/api/app/auto_spatial_advisory/process_hfi.py @@ -14,8 +14,11 @@ from app.db.crud.auto_spatial_advisory import ( save_hfi, get_hfi_classification_threshold, HfiClassificationThresholdEnum, save_run_parameters, get_run_parameters_id) +from app.db.crud.snow import get_last_processed_snow_by_source +from app.db.models.snow import SnowSourceEnum from app.auto_spatial_advisory.classify_hfi import classify_hfi from app.auto_spatial_advisory.run_type import RunType +from app.auto_spatial_advisory.snow import apply_snow_mask from app.geospatial import NAD83_BC_ALBERS from app.auto_spatial_advisory.hfi_pmtiles import get_pmtiles_filepath from app.utils.polygonize import polygonize_in_memory @@ -97,16 +100,19 @@ async def process_hfi(run_type: RunType, run_date: date, run_datetime: datetime, f'for_date:{for_date}' )) return + last_processed_snow = await get_last_processed_snow_by_source(session, SnowSourceEnum.viirs) logger.info('Processing HFI %s for run date: %s, for date: %s', run_type, run_date, for_date) perf_start = perf_counter() - key = get_s3_key(run_type, run_date, for_date) - logger.info(f'Key to HFI in object storage: {key}') + hfi_key = get_s3_key(run_type, run_date, for_date) + logger.info(f'Key to HFI in object storage: {hfi_key}') with tempfile.TemporaryDirectory() as temp_dir: temp_filename = os.path.join(temp_dir, 'classified.tif') - classify_hfi(key, temp_filename) - with polygonize_in_memory(temp_filename, 'hfi', 'hfi') as layer: + classify_hfi(hfi_key, temp_filename) + # Create a snow coverage mask from previously downloaded snow data. + snow_masked_hfi_path = await apply_snow_mask(temp_filename, last_processed_snow[0], temp_dir) + with polygonize_in_memory(snow_masked_hfi_path, 'hfi', 'hfi') as layer: # We need a geojson file to pass to tippecanoe temp_geojson = write_geojson(layer, temp_dir) diff --git a/api/app/auto_spatial_advisory/snow.py b/api/app/auto_spatial_advisory/snow.py new file mode 100644 index 000000000..62030c0f7 --- /dev/null +++ b/api/app/auto_spatial_advisory/snow.py @@ -0,0 +1,117 @@ +import os +import numpy as np +from osgeo import gdal +from app import config +from app.db.models.snow import ProcessedSnow + +SNOW_COVERAGE_WARPED_NAME = 'snow_coverage_warped.tif' +SNOW_COVERAGE_MASK_NAME = 'snow_coverage_mask.tif' +MASKED_HFI_PATH_NAME = 'masked_hfi.tif' + +def classify_snow_mask(snow_path: str, temp_dir: str): + """ + Given a path to snow coverage data, re-classify the data to act as a mask for future HFI processing. + A NDSI (ie. snow coverage) value between 0-100 represent snow coverage. Here we define snow coverage + between 10-100. We need to consult the literature or data scientists on proper use of NDSI. + """ + source = gdal.Open(snow_path, gdal.GA_ReadOnly) + source_band = source.GetRasterBand(1) + source_data = source_band.ReadAsArray() + # In the classified data 0 is assigned to snow covered pixels which will 'cancel' HFI values + # when the rasters are multiplied later on. QA values in the original data are assigned a value + # of 1 so they dont impact HFI calculations for now. + classified = np.where((source_data > 10) & (source_data <= 100), 0, 1) + output_driver = gdal.GetDriverByName("GTiff") + snow_mask_path = os.path.join(temp_dir, SNOW_COVERAGE_MASK_NAME) + snow_mask = output_driver.Create(snow_mask_path, xsize=source_band.XSize, + ysize=source_band.YSize, bands=1, eType=gdal.GDT_Byte) + snow_mask.SetGeoTransform(source.GetGeoTransform()) + snow_mask.SetProjection(source.GetProjection()) + snow_mask_band = snow_mask.GetRasterBand(1) + snow_mask_band.WriteArray(classified) + snow_mask_band = None + snow_mask = None + source_data = None + source_band = None + source = None + return os.path.join(temp_dir, SNOW_COVERAGE_MASK_NAME) + + +async def prepare_snow_mask(hfi_path: str, last_processed_snow: ProcessedSnow, temp_dir: str): + + # Open the HFI tiff to use as a source of parameters for preparation of the snow mask + source = gdal.Open(hfi_path, gdal.GA_ReadOnly) + geo_transform = source.GetGeoTransform() + x_res = geo_transform[1] + y_res = -geo_transform[5] + minx = geo_transform[0] + maxy = geo_transform[3] + maxx = minx + geo_transform[1] * source.RasterXSize + miny = maxy + geo_transform[5] * source.RasterYSize + extent = [minx, miny, maxx, maxy] + source_projection = source.GetProjection() + + bucket = config.get('OBJECT_STORE_BUCKET') + for_date = last_processed_snow.for_date + # The filename of the snow coverage tiff in our object store, prepended with "vsis3" - which tells GDAL to use + # it's S3 virtual file system driver to read the file. + # https://gdal.org/user/virtual_file_systems.html + snow_key = f"/vsis3/{bucket}/snow_coverage/{for_date.strftime('%Y-%m-%d')}/clipped_snow_coverage_{for_date.strftime('%Y-%m-%d')}_epsg4326.tif" + snow_mask_warped_output_path = os.path.join(temp_dir, SNOW_COVERAGE_WARPED_NAME) + # Perform reprojection to Lambert Conformal Conic to match HFI data, crop extent and resample to 2km x 2km pixels + gdal.Warp(snow_mask_warped_output_path, snow_key, dstSRS=source_projection, + outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=gdal.GRA_NearestNeighbour) + source = None + return snow_mask_warped_output_path + +def apply_snow_mask_to_hfi(hfi_path: str, snow_mask_path: str, temp_dir: str): + # Open the hfi data + source_tiff = gdal.Open(hfi_path, gdal.GA_ReadOnly) + source_band = source_tiff.GetRasterBand(1) + source_data = source_band.ReadAsArray() + + # Open the snow mask + snow_tiff = gdal.Open(snow_mask_path, gdal.GA_ReadOnly) + snow_band = snow_tiff.GetRasterBand(1) + snow_data = snow_band.ReadAsArray() + + # The snow mask tiff has values of 0 or 1 where 0 represents areas covered by snow and 1 represents + # snow free areas. Multiply the rasters to apply the mask. + masked_data = np.multiply(source_data, snow_data) + masked_hfi_path = os.path.join(temp_dir, MASKED_HFI_PATH_NAME) + output_driver = gdal.GetDriverByName("GTiff") + # Create an object with the same dimensions as the hfi input + masked_hfi_tiff = output_driver.Create(masked_hfi_path, xsize=source_band.XSize, + ysize=source_band.YSize, bands=1, eType=gdal.GDT_Byte) + # Set the geotransform and projection to the same as the input. + masked_hfi_tiff.SetGeoTransform(source_tiff.GetGeoTransform()) + masked_hfi_tiff.SetProjection(source_tiff.GetProjection()) + + # Write the classified data to the band. + masked_hfi_tiff_band = masked_hfi_tiff.GetRasterBand(1) + masked_hfi_tiff_band.SetNoDataValue(0) + masked_hfi_tiff_band.WriteArray(masked_data) + + # Explicit delete to make sure underlying resources are cleared up! + source_data = None + source_band = None + source_tiff = None + snow_data = None + snow_band = None + snow_tiff = None + masked_hfi_tiff_band = None + masked_hfi_tiff = None + + return masked_hfi_path + + +async def apply_snow_mask(hfi_path: str, last_processed_snow: ProcessedSnow, temp_dir: str): + gdal.SetConfigOption('AWS_SECRET_ACCESS_KEY', config.get('OBJECT_STORE_SECRET')) + gdal.SetConfigOption('AWS_ACCESS_KEY_ID', config.get('OBJECT_STORE_USER_ID')) + gdal.SetConfigOption('AWS_S3_ENDPOINT', config.get('OBJECT_STORE_SERVER')) + gdal.SetConfigOption('AWS_VIRTUAL_HOSTING', 'FALSE') + + snow_mask_path = await prepare_snow_mask(hfi_path, last_processed_snow, temp_dir) + snow_masked_classified_path = classify_snow_mask(snow_mask_path, temp_dir) + hfi_snow_masked_path = apply_snow_mask_to_hfi(hfi_path, snow_masked_classified_path, temp_dir) + return hfi_snow_masked_path