Integrate tile2net task

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

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

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",

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",

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,
+            ),
+        ),
+    ]

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'},
+            }
+        ],
+    },
 }

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
 
 

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()