diff --git a/dev/Dockerfile b/dev/Dockerfile index 6507669a..e2fa8216 100644 --- a/dev/Dockerfile +++ b/dev/Dockerfile @@ -2,11 +2,15 @@ FROM python:3.10-slim # Install system libraries for Python packages: RUN apt-get update && \ apt-get install --no-install-recommends --yes \ - libpq-dev libvips-dev gcc libc6-dev gdal-bin git && \ + libpq-dev libvips-dev gcc libc6-dev gdal-bin git libgl1-mesa-glx && \ rm -rf /var/lib/apt/lists/* -ENV PYTHONDONTWRITEBYTECODE 1 -ENV PYTHONUNBUFFERED 1 +ENV PYTHONDONTWRITEBYTECODE=1 +ENV PYTHONUNBUFFERED=1 + +# Install tile2net (may take up to 30 mins) +RUN git clone https://github.com/VIDA-NYU/tile2net.git +RUN python -m pip install ./tile2net COPY ./setup.py /opt/uvdat-server/setup.py COPY ./manage.py /opt/uvdat-server/manage.py diff --git a/docker-compose.override.yml b/docker-compose.override.yml index 55e4e352..148a01e7 100644 --- a/docker-compose.override.yml +++ b/docker-compose.override.yml @@ -43,6 +43,15 @@ services: env_file: ./dev/.env.docker-compose volumes: - .:/opt/uvdat-server + # NVIDIA driver and increased memory limit are required for tile2net inference + deploy: + resources: + reservations: + devices: + - driver: nvidia + count: 1 + capabilities: [gpu] + shm_size: 1gb depends_on: - postgres - rabbitmq diff --git a/sample_data/use_cases/boston_floods/datasets.json b/sample_data/use_cases/boston_floods/datasets.json index c39782f5..b8baa520 100644 --- a/sample_data/use_cases/boston_floods/datasets.json +++ b/sample_data/use_cases/boston_floods/datasets.json @@ -99,6 +99,19 @@ "trim_distribution_percentage": 0.01 } }, + { + "name": "Boston Orthoimagery", + "description": "Sourced from https://gis.data.mass.gov/maps/6c7009c789354573a42af7251fb768a4/about", + "category": "imagery", + "type": "raster", + "files": [ + { + "url": "https://data.kitware.com/api/v1/item/6744df3c22f196cb5d931161/download", + "path": "boston/boston_orthoimagery.tiff" + } + ], + "style_options": {} + }, { "name": "Boston Neighborhoods", "description": "The Neighborhood boundaries data layer is a combination of zoning neighborhood boundaries, zip code boundaries and 2010 Census tract boundaries", diff --git a/sample_data/use_cases/boston_floods/projects.json b/sample_data/use_cases/boston_floods/projects.json index 409a7911..87a8c874 100644 --- a/sample_data/use_cases/boston_floods/projects.json +++ b/sample_data/use_cases/boston_floods/projects.json @@ -10,6 +10,7 @@ "MBTA Rapid Transit", "MBTA Commuter Rail", "Massachusetts Elevation Data", + "Boston Orthoimagery", "Boston Hurricane Surge Inundation Zones", "Bsoton FEMA National Flood Hazard", "Boston Neighborhoods", diff --git a/uvdat/core/migrations/0008_curbs_aid_simulation.py b/uvdat/core/migrations/0008_curbs_aid_simulation.py new file mode 100644 index 00000000..83eea4b1 --- /dev/null +++ b/uvdat/core/migrations/0008_curbs_aid_simulation.py @@ -0,0 +1,25 @@ +# Generated by Django 5.0.7 on 2024-11-20 22:24 + +from django.db import migrations, models + + +class Migration(migrations.Migration): + + dependencies = [ + ('core', '0007_delete_derivedregion'), + ] + + operations = [ + migrations.AlterField( + model_name='simulationresult', + name='simulation_type', + field=models.CharField( + choices=[ + ('FLOOD_1', 'Flood Scenario 1'), + ('RECOVERY', 'Recovery Scenario'), + ('SEGMENT_CURBS', 'Segment Curbs'), + ], + max_length=13, + ), + ), + ] diff --git a/uvdat/core/models/simulations.py b/uvdat/core/models/simulations.py index 00bd9bdc..6115a23a 100644 --- a/uvdat/core/models/simulations.py +++ b/uvdat/core/models/simulations.py @@ -12,6 +12,7 @@ class SimulationResult(TimeStampedModel): class SimulationType(models.TextChoices): FLOOD_1 = 'FLOOD_1', 'Flood Scenario 1' RECOVERY = 'RECOVERY', 'Recovery Scenario' + SEGMENT_CURBS = 'SEGMENT_CURBS', 'Segment Curbs' simulation_type = models.CharField( max_length=max(len(choice[0]) for choice in SimulationType.choices), @@ -101,4 +102,19 @@ def run(self, **kwargs): }, ], }, + 'SEGMENT_CURBS': { + 'description': """ + Use tile2net to detect roads, sidewalks, footpaths, + and crosswalks in the selected image. + """, + 'output_type': 'dataset', + 'func': uvdat_simulations.segment_curbs, + 'args': [ + { + 'name': 'imagery_layer', + 'type': RasterMapLayer, + 'options_query': {'dataset__category': 'imagery'}, + } + ], + }, } diff --git a/uvdat/core/tasks/map_layers.py b/uvdat/core/tasks/map_layers.py index 31aa342f..00d2d7ab 100644 --- a/uvdat/core/tasks/map_layers.py +++ b/uvdat/core/tasks/map_layers.py @@ -46,10 +46,63 @@ def add_styling(geojson_data, style_options): return geopandas.GeoDataFrame.from_features(features) +def rasterio_convert_raster( + map_layer, temp_dir, raw_data_path, transparency_threshold, trim_distribution_percentage +): + raster_path = Path(temp_dir, 'raster.tiff') + with open(raw_data_path, 'rb') as raw_data: + # read original data with rasterio + input_data = rasterio.open(raw_data) + output_data = rasterio.open( + raster_path, + 'w', + driver='GTiff', + height=input_data.height, + width=input_data.width, + count=1, + dtype=numpy.float32, + crs=input_data.crs, + transform=input_data.transform, + ) + band = input_data.read(1) + + # alter data according to style options + if trim_distribution_percentage: + # trim a number of values from both ends of the distribution + histogram, bin_edges = numpy.histogram(band, bins=1000) + trim_n = band.size * trim_distribution_percentage + new_min = None + new_max = None + sum_values = 0 + for bin_index, bin_count in enumerate(histogram): + bin_edge = bin_edges[bin_index] + sum_values += bin_count + if new_min is None and sum_values > trim_n: + new_min = bin_edge + if new_max is None and sum_values > band.size - trim_n: + new_max = bin_edge + if new_min: + band[band < new_min] = new_min + if new_max: + band[band > new_max] = new_max + if transparency_threshold is not None: + # clamp values below transparency threshold + band[band < transparency_threshold] = transparency_threshold + + band_range = [float(band.min()), float(band.max())] + map_layer.default_style['data_range'] = band_range + + output_data.write(band, 1) + output_data.close() + return raster_path + + def create_raster_map_layer(file_item, style_options): """Save a RasterMapLayer from a FileItem's contents.""" + import large_image import large_image_converter + # create new raster map layer object new_map_layer = RasterMapLayer.objects.create( dataset=file_item.dataset, metadata={}, @@ -58,66 +111,61 @@ def create_raster_map_layer(file_item, style_options): ) print('\t', f'RasterMapLayer {new_map_layer.id} created.') + transparency_threshold = style_options.get('transparency_threshold') + trim_distribution_percentage = style_options.get('trim_distribution_percentage') + with tempfile.TemporaryDirectory() as temp_dir: + raster_path = None raw_data_path = Path(temp_dir, 'raw_data.tiff') + cog_raster_path = Path(temp_dir, 'cog_raster.tiff') + # write original data from file field to file in temp_dir with open(raw_data_path, 'wb') as raw_data: with file_item.file.open('rb') as raw_data_archive: raw_data.write(raw_data_archive.read()) - transparency_threshold = style_options.get('transparency_threshold') - trim_distribution_percentage = style_options.get('trim_distribution_percentage') - - raster_path = Path(temp_dir, 'raster.tiff') - with open(raw_data_path, 'rb') as raw_data: - input_data = rasterio.open(raw_data) - output_data = rasterio.open( - raster_path, - 'w', - driver='GTiff', - height=input_data.height, - width=input_data.width, - count=1, - dtype=numpy.float32, - crs=input_data.crs, - transform=input_data.transform, + if transparency_threshold or trim_distribution_percentage: + # if data needs to be altered according to style options, use rasterio + raster_path = rasterio_convert_raster( + new_map_layer, + temp_dir, + raw_data_path, + transparency_threshold, + trim_distribution_percentage, ) - band = input_data.read(1) - - if trim_distribution_percentage: - # trim a number of values from both ends of the distribution - histogram, bin_edges = numpy.histogram(band, bins=1000) - trim_n = band.size * trim_distribution_percentage - new_min = None - new_max = None - sum_values = 0 - for bin_index, bin_count in enumerate(histogram): - bin_edge = bin_edges[bin_index] - sum_values += bin_count - if new_min is None and sum_values > trim_n: - new_min = bin_edge - if new_max is None and sum_values > band.size - trim_n: - new_max = bin_edge - if new_min: - band[band < new_min] = new_min - if new_max: - band[band > new_max] = new_max - - if transparency_threshold is not None: - band[band < transparency_threshold] = transparency_threshold - - band_range = [float(band.min()), float(band.max())] - new_map_layer.default_style['data_range'] = band_range - - output_data.write(band, 1) - output_data.close() - - cog_raster_path = Path(temp_dir, 'cog_raster.tiff') - large_image_converter.convert(str(raster_path), str(cog_raster_path)) - with open(cog_raster_path, 'rb') as cog_raster_file: - new_map_layer.cloud_optimized_geotiff.save( - cog_raster_path, ContentFile(cog_raster_file.read()) + else: + try: + # if applicable, use large_image to interpret data + source = large_image.open(raw_data_path) + if source.geospatial: + raster_path = raw_data_path + min_val, max_val = None, None + source._scanForMinMax(numpy.int64) + for _frame, range_spec in source._bandRanges.items(): + frame_min = numpy.min(range_spec.get('min', numpy.empty(1))) + frame_max = numpy.max(range_spec.get('max', numpy.empty(1))) + if min_val is None or frame_min < min_val: + min_val = frame_min + if max_val is None or frame_max < max_val: + max_val = frame_max + new_map_layer.default_style['data_range'] = [int(min_val), int(max_val)] + except large_image.exceptions.TileSourceError: + # if original data cannot be interpreted by large_image, use rasterio + raster_path = rasterio_convert_raster( + new_map_layer, + temp_dir, + raw_data_path, + transparency_threshold, + trim_distribution_percentage, ) + # use large_image to convert new raster data to COG + large_image_converter.convert(str(raster_path), str(cog_raster_path)) + with open(cog_raster_path, 'rb') as cog_raster_file: + # save COG to new raster map layer + new_map_layer.cloud_optimized_geotiff.save( + cog_raster_path, ContentFile(cog_raster_file.read()) + ) + return new_map_layer diff --git a/uvdat/core/tasks/simulations.py b/uvdat/core/tasks/simulations.py index d552f8c1..91dcc5f1 100644 --- a/uvdat/core/tasks/simulations.py +++ b/uvdat/core/tasks/simulations.py @@ -1,9 +1,13 @@ +from datetime import datetime +import os from pathlib import Path import random import tempfile +import zipfile from celery import shared_task -from django_large_image import tilesource +from django.core.files.base import ContentFile +from django_large_image import tilesource, utilities import large_image import shapely @@ -117,3 +121,89 @@ def recovery_scenario(simulation_result_id, node_failure_simulation_result, reco 'node_recoveries': node_recoveries, } result.save() + + +@shared_task +def segment_curbs(simulation_result_id, imagery_layer): + from tile2net import Raster + + from uvdat.core.models import Dataset, FileItem, RasterMapLayer, SimulationResult + + result = SimulationResult.objects.get(id=simulation_result_id) + try: + imagery_layer = RasterMapLayer.objects.get(id=imagery_layer) + except Exception: + result.error_message = 'Object not found.' + result.save() + return + + # use django-large-image to open file; + # uses slippy map coords instead of pixel coords to comply with tile2net + imagery_path = utilities.field_file_to_local_path(imagery_layer.cloud_optimized_geotiff) + source = tilesource.get_tilesource_from_path(imagery_path, projection='EPSG:3857') + metadata = source.getMetadata() + zoom = metadata['levels'] - 1 + bounds = metadata.get('sourceBounds') + bbox = [bounds['ymin'], bounds['xmin'], bounds['ymax'], bounds['xmax']] + + with tempfile.TemporaryDirectory() as tmp: + xyz_folder = Path(tmp, 'xyz') + xyz_folder.mkdir(parents=True, exist_ok=True) + output_folder = Path(tmp, 'output') + output_folder.mkdir(parents=True, exist_ok=True) + raster = Raster( + location=bbox, + name='area', + input_dir=f'{str(xyz_folder)}/x_y.png', + output_dir=output_folder, + zoom=zoom, + ) + for tile_row in raster.tiles: + for tile_spec in tile_row: + x, y = tile_spec.xtile, tile_spec.ytile + tile_path = xyz_folder / f'{x}_{y}.png' + tile = source.getTile(x, y, zoom, pilImageAllowed=True) + tile.save(tile_path) + + try: + stitch_step = 4 + raster.generate(stitch_step) + # inference step requires NVIDIA driver and increased memory limit + raster.inference() + + dataset_ids = [] + for result_set in ['polygons', 'network']: + result_folder = next(output_folder.glob(f'area/{result_set}')) + zip_path = output_folder / f'{result_set}.zip' + if result_folder.exists(): + with zipfile.ZipFile(zip_path, 'w', zipfile.ZIP_DEFLATED) as z: + for result_file_path in result_folder.glob('**/*'): + if result_file_path.is_file(): + z.write(str(result_file_path)) + if zip_path.exists(): + dataset_name = f'Generated {result_set} for {imagery_layer.dataset.name}' + existing_count = Dataset.objects.filter(name__contains=dataset_name).count() + if existing_count: + dataset_name += f' ({existing_count + 1})' + dataset = Dataset.objects.create( + name=dataset_name, + description='Segmentation generated by tile2net from orthoimagery', + category='segmentation', + dataset_type=Dataset.DatasetType.VECTOR, + metadata={'creation_time': datetime.now().strftime('%Y-%m-%d %H:%M:%S')}, + ) + result.project.datasets.add(dataset) + file_item = FileItem.objects.create( + name=zip_path.name, + dataset=dataset, + file_type='zip', + file_size=os.path.getsize(zip_path), + ) + with zip_path.open('rb') as f: + file_item.file.save(zip_path, ContentFile(f.read())) + dataset.spawn_conversion_task(asynchronous=False) + dataset_ids.append(dataset.id) + result.output_data = {'dataset_ids': dataset_ids} + except Exception as e: + result.error_message = str(e) + result.save() diff --git a/web/src/components/DatasetList.vue b/web/src/components/DatasetList.vue index 494ea42d..7ad0df00 100644 --- a/web/src/components/DatasetList.vue +++ b/web/src/components/DatasetList.vue @@ -225,7 +225,7 @@ onMounted(fetchAllLayers); :model-value="selectedIds.includes(dataset.id)" @update:model-value="toggleDatasets($event, [dataset])" /> - {{ dataset.name }} +
{{ dataset.name }}
diff --git a/web/src/components/OptionsDrawerContents.vue b/web/src/components/OptionsDrawerContents.vue index 8a34578f..0a1b4211 100644 --- a/web/src/components/OptionsDrawerContents.vue +++ b/web/src/components/OptionsDrawerContents.vue @@ -50,7 +50,7 @@ export default { VectorDatasetLayer | RasterDatasetLayer | undefined >(); const opacity = ref(1); - const colormap = ref("terrain"); + const colormap = ref(); const layerRange = ref([]); const colormapRange = ref([]); const applyToAll = ref(true); @@ -92,7 +92,7 @@ export default { function populateRefs() { if (!currentDatasetLayer.value) { opacity.value = 1; - colormap.value = "terrain"; + colormap.value = undefined; layerRange.value = []; colormapRange.value = []; @@ -129,7 +129,7 @@ export default { } // Populate raster refs, if necessary - colormap.value = defaultStyle.palette || "terrain"; + colormap.value = defaultStyle.palette || undefined; layerRange.value = defaultStyle.data_range?.map((v: number) => Math.round(v)) || []; colormapRange.value = layerRange.value; @@ -194,12 +194,15 @@ export default { return; } - styleRasterDatasetLayer(currentDatasetLayer.value, { - colormap: { - palette: colormap.value, - range: colormapRange.value, - }, - }); + const style = colormap.value + ? { + colormap: { + palette: colormap.value, + range: colormapRange.value, + }, + } + : {}; + styleRasterDatasetLayer(currentDatasetLayer.value, style); } function getNetworkNodeName(nodeId: number) { @@ -309,6 +312,7 @@ export default { diff --git a/web/src/components/SimulationsPanel.vue b/web/src/components/SimulationsPanel.vue index cd39ae75..ff4abac0 100644 --- a/web/src/components/SimulationsPanel.vue +++ b/web/src/components/SimulationsPanel.vue @@ -4,8 +4,14 @@ import { currentSimulationType, currentProject, selectedDatasetLayers, + availableDatasets, } from "@/store"; -import { getSimulationResults, runSimulation } from "@/api/rest"; +import { + getDatasetLayers, + getProjectDatasets, + getSimulationResults, + runSimulation, +} from "@/api/rest"; import NodeAnimation from "./NodeAnimation.vue"; import { isRasterDatasetLayer, @@ -119,6 +125,23 @@ export default { }[]; } + async function populateActiveResultOutputs() { + if (activeResult.value && currentProject.value) { + availableDatasets.value = await getProjectDatasets( + currentProject.value.id + ); + const datasetIds = activeResult.value.output_data?.dataset_ids; + activeResult.value.output_data.datasets = await Promise.all( + availableDatasets.value + .filter((d) => datasetIds.includes(d.id)) + .map(async (d) => { + d.map_layers = await getDatasetLayers(d.id); + return d; + }) + ); + } + } + function pollForActiveResultOutput() { if (!availableResults.value) { clearInterval(outputPoll.value); @@ -154,6 +177,7 @@ export default { watch(activeResult, () => { if (activeResult.value) { populateActiveResultInputs(); + populateActiveResultOutputs(); if (!outputPoll.value && !activeResult.value.output_data) { outputPoll.value = setInterval(pollForActiveResultOutput, 3000); } @@ -288,6 +312,29 @@ export default { :nodeRecoveries="result.output_data.node_recoveries" /> +
+ + + + {{ dataset.name }} + + + Show on Map + + + + + +
Unknown simulation output type {{ currentSimulationType.output_type }} diff --git a/web/src/layers.ts b/web/src/layers.ts index 2b7c71d3..d55f2520 100644 --- a/web/src/layers.ts +++ b/web/src/layers.ts @@ -505,17 +505,16 @@ export function styleRasterDatasetLayer( | Record | undefined; const colormapPalette = (options?.colormap?.palette || - defaultStyle?.palette || - "terrain") as string; + defaultStyle?.palette) as string; const colormapRange = options?.colormap?.range || defaultStyle?.data_range; const nodataValue = options.nodata || defaultStyle?.transparency_threshold; const tileParams: Record = { projection: "EPSG:3857", - band: "1", }; if (colormapPalette) { tileParams.palette = colormapPalette; + tileParams.band = "1"; } if (Array.isArray(colormapRange) && colormapRange?.length === 2) { tileParams.min = colormapRange[0]; diff --git a/web/src/types.ts b/web/src/types.ts index 2326145c..33ad5715 100644 --- a/web/src/types.ts +++ b/web/src/types.ts @@ -268,6 +268,8 @@ export interface SimulationResult { output_data: { node_failures: []; node_recoveries: []; + datasets: Dataset[]; + dataset_ids: number[]; }; error_message: string; created: string;