From 96494fe0b199b46744add4bbf66a8a34f01d7af7 Mon Sep 17 00:00:00 2001 From: Terence Tuhinanshu Date: Fri, 28 Jun 2024 17:34:03 -0400 Subject: [PATCH 1/6] Add dependencies for fetching Global IO LULC data --- src/mmw/requirements/base.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/mmw/requirements/base.txt b/src/mmw/requirements/base.txt index a44baff8f..b7c38d35c 100644 --- a/src/mmw/requirements/base.txt +++ b/src/mmw/requirements/base.txt @@ -30,3 +30,6 @@ redis<4 numpy==1.22.* dictdiffer==0.9.0 BMPxlsx==3.0.0 +rasterio==1.3.10 +pystac-client==0.8.2 +shapely==2.0.4 From c2c942c16038a3ef41f9dfa69f08aed116fd8e2d Mon Sep 17 00:00:00 2001 From: Terence Tuhinanshu Date: Thu, 11 Jul 2024 17:03:41 -0400 Subject: [PATCH 2/6] Add Celery task to query STAC and generate histogram This is based on preparatory work done in https://github.com/rajadain/mmw-io-10m-lulc-summary/blob/main/main.py We create a Celery task that takes a STAC URL, collection, asset, and optional filters, and generates a histogram from it. There is an additional method to format this into an MMW Geoprocessing like shape, which allows for compatibility and modularity. --- src/mmw/apps/geoprocessing_api/stac.py | 179 +++++++++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 src/mmw/apps/geoprocessing_api/stac.py diff --git a/src/mmw/apps/geoprocessing_api/stac.py b/src/mmw/apps/geoprocessing_api/stac.py new file mode 100644 index 000000000..fe843c435 --- /dev/null +++ b/src/mmw/apps/geoprocessing_api/stac.py @@ -0,0 +1,179 @@ +from celery import shared_task +from pystac_client import Client +from rasterio.coords import BoundingBox +from rasterio.enums import Resampling +from rasterio.mask import mask +from rasterio.merge import merge +from rasterio.warp import ( + calculate_default_transform, + transform_geom, + reproject, +) +from shapely.geometry import box, mapping, shape +import json +import numpy as np +import os +import rasterio as rio +import tempfile + + +@shared_task +def query_histogram(geojson, url, collection, asset, filter): + # TODO Validate inputs + + aoi = shape(json.loads(geojson)) + + # Get list of intersecting tiffs from catalog + client = Client.open(url) + search = client.search( + collections=collection, + intersects=mapping(aoi), + filter=filter, + ) + tiffs = [item.assets[asset].href for item in search.items()] + + # Sort the tiffs, since the merging is order sensitive. We use the "first" + # method, using the first image for overlapping pixels. + tiffs = sorted(tiffs) + + # TODO Handle empty tiff list + + # Find the Albers Equal Area CRS for this AoI + dst_crs = get_albers_crs_for_aoi(aoi) + + # Reproject the tiffs and clip to the AoI to make tiles + clips = [] + for tiff in tiffs: + clip_data, clip_transform, clip_meta = clip_and_reproject_tile( + tiff, aoi, dst_crs + ) + temp_file = tempfile.NamedTemporaryFile(delete=False, suffix='.tif') + with rio.open(temp_file.name, 'w', **clip_meta) as dst: + dst.write(clip_data) + clips.append(temp_file.name) + + # Merge the clipped rasters + datasets = [rio.open(clip_path) for clip_path in clips] + merged_data, merged_transform = merge(datasets, method='first') + + # Count how many of each class type there are + values, counts = np.unique(merged_data, return_counts=True) + histogram = list(zip(values.tolist(), counts.tolist())) + + # Ensure the result is JSON serializable + histogram = {int(value): int(count) for value, count in histogram} + + # Calculate pixel size from dataset resolution + pixel_size = datasets[0].res[0] * datasets[0].res[1] + + # Close datasets + for ds in datasets: + ds.close() + + # Clean up temporary files + for temp_file in clips: + os.remove(temp_file) + + return { + 'result': histogram, + 'pixel_size': pixel_size, + } + + +@shared_task +def format_as_mmw_geoprocessing(result): + """ + Return results formatted as MMW Geoprocessing + + To play well with existing patterns, and to allow for switching this out + for other implementations in the future. + """ + # Pass errors on + if 'error' in result: + return result + + formatted_histogram = { + f'List({v})': int(c) for v, c in result['result'].items() + } + result['result'] = formatted_histogram + + return result + + +def get_albers_crs_for_aoi(aoi): + """ + Return the Albers Equal Area Projection for the given AoI + + Since we want to calculate area, we need to use an Equal Area projection, + but this differs based on where you are in the globe. We use rough bounding + boxes to see if the shape neatly fits within one of the continents. If not, + we fall back to a global approximation. + """ + + if aoi.within(box(-170, 15, -50, 75)): # North America + return 'EPSG:5070' + elif aoi.within(box(-10, 34, 40, 72)): # Europe + return 'EPSG:3035' + elif aoi.within(box(25, -10, 180, 60)): # Asia + return 'EPSG:102025' + elif aoi.within(box(-20, -35, 55, 38)): # Africa + return 'EPSG:102022' + elif aoi.within(box(-90, -60, -30, 15)): # South America + return 'EPSG:102033' + elif aoi.within(box(112, -45, 155, -10)): # Australia + return 'EPSG:102034' + else: # Global + return 'EPSG:54017' # Behrmann Equal Area Cylindrical + + +def clip_and_reproject_tile(tiff, aoi, dst_crs, resampling=Resampling.nearest): + """ + Clip a tiff to the given AoI and reproject to destination CRS + + We clip the tiff first, and then reproject only the clipped bits to + minimize transformation. Uses nearest neighbor resampling by default. + + Returns the data, transform, and meta after clipping and reprojecting. + """ + aoi_crs = aoi.crs if hasattr(aoi, 'crs') else 4326 + + with rio.open(tiff) as src: + # Reproject the AoI into the tiff's CRS, then clip the tiff to + # the reprojected AoI's bounding box + reprojected_aoi = transform_geom(aoi_crs, src.crs, aoi) + reprojected_aoi_bbox = BoundingBox(*shape(reprojected_aoi).bounds) + clip_data, clip_transform = mask(src, [reprojected_aoi], crop=True) + + # Define the output metadata for the reprojected clip + dst_transform, width, height = calculate_default_transform( + src.crs, + dst_crs, + clip_data.shape[2], + clip_data.shape[1], + *reprojected_aoi_bbox, + ) + dst_meta = src.meta.copy() + dst_meta.update( + { + 'crs': dst_crs, + 'transform': dst_transform, + 'width': width, + 'height': height, + } + ) + + # Reproject the clipped data to the destination CRS + reprojected_clip_data = np.empty( + (src.count, height, width), dtype=src.dtypes[0] + ) + reproject( + source=clip_data, + destination=reprojected_clip_data, + src_transform=clip_transform, + src_crs=src.crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=resampling, + ) + + return reprojected_clip_data, dst_transform, dst_meta From 972c74ca3cb7bd7d77c38054938530860f6e2ac2 Mon Sep 17 00:00:00 2001 From: Terence Tuhinanshu Date: Fri, 12 Jul 2024 11:56:20 -0400 Subject: [PATCH 3/6] Add IO LULC Class Map --- src/mmw/mmw/settings/layer_classmaps.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/mmw/mmw/settings/layer_classmaps.py b/src/mmw/mmw/settings/layer_classmaps.py index c4bad6ec9..f315f7041 100644 --- a/src/mmw/mmw/settings/layer_classmaps.py +++ b/src/mmw/mmw/settings/layer_classmaps.py @@ -51,3 +51,18 @@ 11: ['con_ease', 'Conservation Easement'], 12: ['ag_ease', 'Agricultural Easement'], } + +# Map for the Impact Observatory Land Use Land Classification scheme +# https://www.impactobservatory.com/maps-for-good/ +IO_LULC = { + 0: ['no_data', 'No Data'], + 1: ['water', 'Water'], + 2: ['trees', 'Trees'], + 4: ['flooded_vegetation', 'Flooded vegetation'], + 5: ['crops', 'Crops'], + 7: ['built_area', 'Built area'], + 8: ['bare_ground', 'Bare ground'], + 9: ['snow_ice', 'Snow/ice'], + 10: ['clouds', 'Clouds'], + 11: ['rangeland', 'Rangeland'], +} From ebbf3b38ded6554c4becd67f7bd5a019e065b341 Mon Sep 17 00:00:00 2001 From: Terence Tuhinanshu Date: Mon, 15 Jul 2024 10:26:04 -0400 Subject: [PATCH 4/6] Add analyze task for global land cover This task takes the output of the STAC Histogram Query and calculates the area and coverage percent from it, and formats it into the standard MMW output. We remove IO Class 0 NOADAT from the standard set of outputs so that we don't include it in the final output. --- src/mmw/apps/geoprocessing_api/tasks.py | 40 +++++++++++++++++++++++++ src/mmw/mmw/settings/layer_classmaps.py | 2 +- 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/src/mmw/apps/geoprocessing_api/tasks.py b/src/mmw/apps/geoprocessing_api/tasks.py index 6377d2055..e7e2c9a1f 100644 --- a/src/mmw/apps/geoprocessing_api/tasks.py +++ b/src/mmw/apps/geoprocessing_api/tasks.py @@ -184,6 +184,46 @@ def area(dictionary, key, default=0): } +@shared_task(throws=Exception) +def analyze_global_land(result, year): + if 'error' in result: + raise Exception(f'[analyze_global_land_{year}] {result["error"]}') + + histogram = parse(result['result']) + pixel_size = result['pixel_size'] + + # Count up all the legitimate IO LULC classes + # This excludes NODATA which is IO LULC 0 + total_count = sum( + [ + count + for ioclass, count in histogram.items() + if ioclass in layer_classmaps.IO_LULC + ] + ) + + categories = [] + + for ioclass, (code, name) in layer_classmaps.IO_LULC.items(): + categories.append( + { + 'area': histogram.get(ioclass, 0) * pixel_size, + 'code': code, + 'coverage': float(histogram.get(ioclass, 0)) / total_count, + 'ioclass': ioclass, + 'type': name, + } + ) + + return { + 'survey': { + 'name': f'global_land_io_{year}', + 'displayName': f'Global Land Use/Cover {year}', + 'categories': categories, + } + } + + @shared_task(throws=Exception) def analyze_soil(result, area_of_interest=None): if 'error' in result: diff --git a/src/mmw/mmw/settings/layer_classmaps.py b/src/mmw/mmw/settings/layer_classmaps.py index f315f7041..7523066a8 100644 --- a/src/mmw/mmw/settings/layer_classmaps.py +++ b/src/mmw/mmw/settings/layer_classmaps.py @@ -55,7 +55,7 @@ # Map for the Impact Observatory Land Use Land Classification scheme # https://www.impactobservatory.com/maps-for-good/ IO_LULC = { - 0: ['no_data', 'No Data'], + # 0: ['no_data', 'No Data'], 1: ['water', 'Water'], 2: ['trees', 'Trees'], 4: ['flooded_vegetation', 'Flooded vegetation'], From 3632c34178e6b87895fda2e2abdb9935fcf5c493 Mon Sep 17 00:00:00 2001 From: Terence Tuhinanshu Date: Mon, 15 Jul 2024 10:32:13 -0400 Subject: [PATCH 5/6] Add view for analyzing global land This follows the format of the Analyze Land endpoint. Input validation and caching is deferred to a later commit. --- src/mmw/apps/geoprocessing_api/urls.py | 3 + src/mmw/apps/geoprocessing_api/views.py | 132 +++++++++++++++++++++++- 2 files changed, 134 insertions(+), 1 deletion(-) diff --git a/src/mmw/apps/geoprocessing_api/urls.py b/src/mmw/apps/geoprocessing_api/urls.py index a7957ee3a..6b518c77d 100644 --- a/src/mmw/apps/geoprocessing_api/urls.py +++ b/src/mmw/apps/geoprocessing_api/urls.py @@ -12,6 +12,9 @@ name="authtoken"), re_path(r'analyze/land/(?P\w+)/?$', views.start_analyze_land, name='start_analyze_land'), + re_path(r'analyze/global-land/(?P\w+)/?$', + views.start_analyze_global_land, + name='start_analyze_global_land'), re_path(r'analyze/soil/$', views.start_analyze_soil, name='start_analyze_soil'), re_path(r'analyze/animals/$', views.start_analyze_animals, diff --git a/src/mmw/apps/geoprocessing_api/views.py b/src/mmw/apps/geoprocessing_api/views.py index 66f8a5d17..c92c6babc 100644 --- a/src/mmw/apps/geoprocessing_api/views.py +++ b/src/mmw/apps/geoprocessing_api/views.py @@ -37,7 +37,7 @@ from apps.modeling.serializers import AoiSerializer from apps.modeling.views import _initiate_subbasin_gwlfe_job_chain -from apps.geoprocessing_api import exceptions, schemas, tasks +from apps.geoprocessing_api import exceptions, schemas, stac, tasks from apps.geoprocessing_api.calcs import huc12s_for_huc, huc12_for_point from apps.geoprocessing_api.permissions import AuthTokenSerializerAuthentication # noqa from apps.geoprocessing_api.throttling import (BurstRateThrottle, @@ -443,6 +443,136 @@ def start_analyze_land(request, nlcd_year, format=None): ], area_of_interest, user, messages=msg) +@swagger_auto_schema(method='post', + request_body=schemas.ANALYZE_REQUEST, + responses={200: schemas.JOB_STARTED_RESPONSE}) +@decorators.api_view(['POST']) +@decorators.authentication_classes((SessionAuthentication, + TokenAuthentication, )) +@decorators.permission_classes((IsAuthenticated, )) +@decorators.throttle_classes([BurstRateThrottle, SustainedRateThrottle]) +@log_request +def start_analyze_global_land(request, year, format=None): + """ + Starts a job to produce a land-use histogram for a given area and year. + + This endpoint supports year values of 2017, 2018, 2019, 2020, 2021, 2022, + and 2023. + + This endpoint supports global data, as opposed to analyze_land which only + works with the Continental United States (CONUS). The underlying datasource + is Impact Observatory's Global Annual LULC layer hosted on AWS Open Registry. # NOQA + + For more information, see the [technical documentation](https://wikiwatershed.org/documentation/mmw-tech/#overlays-tab-coverage). # NOQA + + ## Response + + You can use the URL provided in the response's `Location` + header to poll for the job's results. + + + **Example of a completed job's `result`** + + +
+ + { + "survey": { + "displayName": "Global Land Use/Cover 2019", + "name": "global_land_io_2019", + "categories": [ + { + "area": 14941947.428126818, + "code": "water", + "coverage": 0.02144532388706757, + "ioclass": 1, + "type": "Water" + }, + { + "area": 109588951.93420613, + "code": "trees", + "coverage": 0.15728676465889277, + "ioclass": 2, + "type": "Trees" + }, + { + "area": 48464.99853478383, + "code": "flooded_vegetation", + "coverage": 0.00006955904481421345, + "ioclass": 4, + "type": "Flooded vegetation" + }, + { + "area": 29756139.065063003, + "code": "crops", + "coverage": 0.042707287182504744, + "ioclass": 5, + "type": "Crops" + }, + { + "area": 488978611.8600828, + "code": "built_area", + "coverage": 0.7018030785899226, + "ioclass": 7, + "type": "Built area" + }, + { + "area": 535341.2912358101, + "code": "bare_ground", + "coverage": 0.0007683447847676015, + "ioclass": 8, + "type": "Bare ground" + }, + { + "area": 0, + "code": "snow_ice", + "coverage": 0, + "ioclass": 9, + "type": "Snow/ice" + }, + { + "area": 0, + "code": "clouds", + "coverage": 0, + "ioclass": 10, + "type": "Clouds" + }, + { + "area": 52896720.20292212, + "code": "rangeland", + "coverage": 0.07591964185203046, + "ioclass": 11, + "type": "Rangeland" + } + ] + } + } + +
+ """ + user = request.user if request.user.is_authenticated else None + area_of_interest, wkaoi, msg = _parse_analyze_input(request) + + filter = { + "op": "like", + "args": [{"property": "id"}, f"%-{year}"], + } + + # TODO Implement caching + + return start_celery_job([ + stac.query_histogram.s( + area_of_interest, + url='https://api.impactobservatory.com/stac-aws', + collection='io-10m-annual-lulc', + asset='supercell', + filter=filter + ), + stac.format_as_mmw_geoprocessing.s(), + tasks.analyze_global_land.s(year), + ], area_of_interest, user, messages=msg) + + @swagger_auto_schema(method='post', request_body=schemas.ANALYZE_REQUEST, responses={200: schemas.JOB_STARTED_RESPONSE}) From fffe8cd4400ea5fe50ed8f163d73d1cb65ff467d Mon Sep 17 00:00:00 2001 From: Terence Tuhinanshu Date: Tue, 16 Jul 2024 10:03:03 -0400 Subject: [PATCH 6/6] Add legend for IO Global LULC land cover --- src/mmw/mmw/settings/layer_settings.py | 6 +++--- src/mmw/sass/utils/_variables.scss | 18 ++++++++++++++++++ 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/src/mmw/mmw/settings/layer_settings.py b/src/mmw/mmw/settings/layer_settings.py index 4392ec4ec..af5a01b62 100644 --- a/src/mmw/mmw/settings/layer_settings.py +++ b/src/mmw/mmw/settings/layer_settings.py @@ -16,7 +16,7 @@ import json from django.contrib.gis.geos import GEOSGeometry -from mmw.settings.layer_classmaps import NLCD, SOIL, PROTECTED_LANDS +from mmw.settings.layer_classmaps import NLCD, SOIL, PROTECTED_LANDS, IO_LULC # [01, 02, ...] style list for layer time sliders MONTH_CODES = [str(m).zfill(2) for m in range(1, 13)] @@ -98,7 +98,7 @@ { 'display': 'IO Global LULC 2023', 'code': 'io-lulc-2023', - 'css_class_prefix': 'io-lulc-2023 nlcd', + 'css_class_prefix': 'io-lulc-2023 io-lulc', 'short_display': 'IO LULC 2023', 'helptext': 'Global land use/land cover dataset produced by ' 'Impact Observatory, Microsoft, and Esri, derived from ' @@ -108,7 +108,7 @@ 'maxZoom': 18, 'opacity': 0.618, 'has_opacity_slider': True, - 'legend_mapping': { key: names[1] for key, names in NLCD.items()}, + 'legend_mapping': { key: names[1] for key, names in IO_LULC.items() }, 'big_cz': True, }, { diff --git a/src/mmw/sass/utils/_variables.scss b/src/mmw/sass/utils/_variables.scss index 9361cdb58..d18c8d79f 100644 --- a/src/mmw/sass/utils/_variables.scss +++ b/src/mmw/sass/utils/_variables.scss @@ -78,6 +78,24 @@ $nlcdColors: ( } } +$ioLulcColors: ( + 1: #419bdf, // Water + 2: #397d49, // Trees + 4: #7a87c6, // Flooded vegetation + 5: #e49635, // Crops + 7: #c4281b, // Built area + 8: #a59b8f, // Bare ground + 9: #a8ebff, // Snow/ice + 10: #616161, // Clouds + 11: #e3e2c3, // Rangeland +); + +@each $id, $color in $ioLulcColors { + .io-lulc-#{$id} { + background-color: $color; + } +} + $protectedLandsColors: ( pra_f: #AFB42B, /* Park or Recreation Area - Federal */ pra_s: #C0CA33, /* Park or Recreation Area - State */