Notebook

In this notebook, you will be guided to :

Create a project and upload an orthomosaic and DSM rasters
Generate a tile layer URLs
Visualize an orthomosaic and a DSM in the notebook
Publish it on arcgisonline

You may execute this notebook one cell after the other with `Shift + Enter`¶

For more options, explore Cell in the top menu bar, or read a tutorial

Requirements¶

delair.ai Python SDK
folium
arcgis
shapely
pyproj

⚠️ Restart the Jupyter kernel after first install

In [ ]:

!pip install python-delairstack folium arcgis shapely pyproj

In [ ]:

import os
from delairstack import DelairStackSDK

In [ ]:

import getpass
platform_url = 'https://www.delair.ai'
login = input('Enter your email ')
password = getpass.getpass('Enter your password ')

In [ ]:

sdk = DelairStackSDK(url=platform_url, user=login, password=password)

Download sample files (images, mesh, raster...)¶

An archive Saint-Papoul.zip containing sample files will be downloaded (if not found in the current directory).

In [ ]:

import urllib.request
import zipfile

try: working_dir
except NameError: working_dir = os.getcwd()
%cd {working_dir}

if not os.path.exists('Saint-Papoul'):
    print('"Saint-Papoul" folder not found')
    if not os.path.exists('Saint-Papoul.zip'):
        print('"Saint-Papoul.zip" not found')
        print('Downloading it...', end=' ')
        url = 'https://delair-transfer.s3-eu-west-1.amazonaws.com/sdks/sample-data/Saint-Papoul.zip'
        filename, _ = urllib.request.urlretrieve(url, 'Saint-Papoul.zip')
        print('OK')

    print('Extracting "Saint-Papoul.zip"...', end=' ')
    with zipfile.ZipFile(filename, 'r') as zip_ref:
        zip_ref.extractall('./Saint-Papoul')
    print('OK')
else:
    print('"Saint-Papoul" folder found. No need to download it again.')

In [ ]:

sample_path = './Saint-Papoul'

In [ ]:

%cd {sample_path}
!ls .

Create the project¶

In [ ]:

my_project = sdk.projects.create(
    name='SDK Map Tiles Tutorial',
    geometry={"coordinates": [[
                [2.0362935018049213, 43.33793077500704],
                [2.052168442639523, 43.33793077500704],
                [2.052168442639523, 43.34757088135606],
                [2.0362935018049213, 43.34757088135606],
                [2.0362935018049213, 43.33793077500704]]],
              "type": "Polygon"})

In [ ]:

print('We just created the project {!r} with id {!r}'.format(
    my_project.name, my_project.id))

Create and upload the orthomosaic¶

In [ ]:

ortho_dataset = sdk.datasets.create_raster_dataset(
  name='Orthomosaic',
  project=my_project.id,
  dataset_format='geotiff',
  categories=['orthomosaic'])

sdk.datasets.upload_file(
  dataset=ortho_dataset.id,
  component='raster',
  file_path='Orthomosaic.tif')

DSM (Digital Surface Model)¶

In [ ]:

dsm_dataset = sdk.datasets.create_raster_dataset(
  name='DSM',
  project=my_project.id,
  dataset_format='geotiff',
  categories=['dsm'])

sdk.datasets.upload_file(
  dataset=dsm_dataset.id,
  component='raster',
  file_path='DSM.tif')

After the raster uploads, they need to be tiled. It can take up to 5 minutes.¶

In [ ]:

from time import sleep
print('Please wait till the rasters are tiled (it can takes up to 5 minutes)...')
while True:
    ds_list = sdk.datasets.describe([ortho_dataset.id, dsm_dataset.id])
    ingestion_statuses = [ds.ingestion.get('status') for ds in ds_list]
    if ingestion_statuses == ['completed', 'completed']:
        break
    else:
        sleep(10)

print('OK 👍')
print('Orthomosaic and DSM have been ingested properly')

Tile layer URLs Generation¶

delair.ai provides a tile service for all raster datasets.

It follows the TMS standard (for example: https://delair.ai/tileserver/tiles/DATASET_ID/%7Bz%7D/%7Bx%7D/%7By%7D.png).

This service is protected using a delair.ai token.

The Python SDK provides the sdk.datasets.share_tiles function:

to quickly generate a token authorized for a specific dataset
to build a URL for the tile service of this dataset.

In [ ]:

ortho_tile_url = sdk.datasets.share_tiles(dataset=ortho_dataset.id)
ortho_tile_url

In [ ]:

dsm_tile_url = sdk.datasets.share_tiles(dataset=dsm_dataset.id)
dsm_tile_url

Determine the center of the raster and its bounding box¶

In [ ]:

# Request the updated dataset properties (such as its geometry)
ortho_dataset = sdk.datasets.describe(ortho_dataset.id)

In [ ]:

from shapely.geometry import shape

def get_center_coords(dataset):
    shapely_geometry = shape(dataset.geometry)
    return list(shapely_geometry.centroid.coords)[0]

center_longitude, center_latitude = get_center_coords(ortho_dataset)
map_center = (center_latitude, center_longitude)
print('Coordinates of the center: lat={:.4f}, long={:.4f}'.format(center_latitude, center_longitude))

In [ ]:

def get_bbox(dataset):
    shapely_geometry = shape(dataset.geometry)
    return shapely_geometry.bounds  # (minx, miny, maxx, maxy)

bbox = get_bbox(ortho_dataset)
print('Bounding box: {}'.format(bbox))

Visualization with Folium¶

Orthomosaic¶

In [ ]:

import folium

m = folium.Map(location=map_center, zoom_start=18)

folium.raster_layers.TileLayer(
    tiles=ortho_tile_url,
    attr='delair.ai',
    max_zoom=20,
    overlay=False,
    control=True,
    bounds=[[bbox[1], bbox[0]], [bbox[3], bbox[2]]]
).add_to(m)
m

In [ ]:

dsm_map = folium.Map(location=map_center, zoom_start=18)

folium.raster_layers.TileLayer(
    tiles=dsm_tile_url,
    attr='delair.ai',
    max_zoom=20,
    overlay=False,
    control=True,
    bounds=[[bbox[1], bbox[0]], [bbox[3], bbox[2]]]
).add_to(dsm_map)
dsm_map

In [ ]:

# Functions to generate a custom tile URL and render it

from urllib.parse import urlencode, urlsplit, parse_qs, urlunsplit
from ipywidgets import interact, fixed, FloatSlider, IntSlider

def add_query_params_to_url(url, param_dict):
    scheme, netloc, path, query_string, fragment = urlsplit(url)
    params = parse_qs(query_string)
    params.update(param_dict)
    new_query_string = urlencode(params, doseq=True)
    return urlunsplit((scheme, netloc, path, new_query_string, fragment))

def generate_custom_tile_url(original_tile_url, **kwargs):
    params = kwargs
    tile_url = add_query_params_to_url(original_tile_url, params)
    return tile_url

def render_folium_map(original_tile_url, **kwargs):
    tile_url = generate_custom_tile_url(original_tile_url, **kwargs)
    m = folium.Map(location=map_center, zoom_start=18)
    folium.raster_layers.TileLayer(
        tiles=tile_url,
        attr='delair.ai',
        max_zoom=20,
        overlay=False,
        control=True,
        bounds=[[bbox[1], bbox[0]], [bbox[3], bbox[2]]]
    ).add_to(m)
    return m

def interactive_dsm_map():
    original_tile_url = dsm_tile_url
    ds = sdk.datasets.describe(dsm_dataset.id)
    min_value = ds.bands[0].get('stats').get('min')
    max_value = ds.bands[0].get('stats').get('max')

    interact(
        render_folium_map,
        original_tile_url=fixed(original_tile_url),
        band_id=fixed(1),
        colormap=['spectral', 'RdYlGn', 'GnYlRd', 'greyscale'],
        outofrange_type=['clip', 'transparent'],
        slope=IntSlider(value=100, max=100, min=0, step=5, continuous_update=False),
        min=FloatSlider(min_value, min=min_value, max=max_value, step=1e-1, continuous_update=False),
        max=FloatSlider(max_value, min=min_value, max=max_value, step=1e-1, continuous_update=False)
    )

For exemple, you can change :

slope = 50
colormap = GnYIRd
min_value around 188
outofrange_type = transparent

In [ ]:

interactive_dsm_map()

With Argis¶

In [ ]:

default_arcgis_url = ''
arcgis_url = input('Enter your arcgis url (for example "https://mycompany.maps.arcgis.com")')
arcgis_username = input('Enter your arcgis username ')
arcgis_password = getpass.getpass('Enter your arcgis password ')

In [ ]:

from arcgis.gis import GIS

print('Connecting to {!r} with the supplied credentials...'.format(arcgis_url))
gis = GIS(arcgis_url, arcgis_username, arcgis_password)
print('OK 👍')

New map¶

You may refer to the ArcGIS documentation and the ArcGIS python package API to get a deeper understanding of the following concepts.

In [ ]:

from pyproj import Proj, transform

def convert_coords_to_mercator(long, lat):
    origin_projection = Proj('epsg:4326')
    mercator_projection = Proj('epsg:3857')
    return transform(origin_projection, mercator_projection, long, lat)

def get_mercator_full_extent(dataset):
    xmin_4326, ymin_4326, xmax_4326, ymax_4326 = get_bbox(dataset)
    xmin, ymin = convert_coords_to_mercator(ymin_4326, xmin_4326)
    xmax, ymax = convert_coords_to_mercator(ymax_4326, xmax_4326)
    
    full_extent = {
        'xmin': xmin,
        'ymin': ymin,
        'xmax': xmax,
        'ymax': ymax,
        'spatialReference': {'wkid': 102100}
    }
    return full_extent

def create_arcgis_layer(name, url):
    proper_url = url.replace('{z}', '{level}').replace('{x}', '{col}').replace('{y}', '{row}')
    return {
        "opacity": 1,
        "visibility": True,
        "title": name,
        "type": "WebTiledLayer",
        "layerType": "WebTiledLayer",
        "templateUrl": proper_url,
        "copyright": "Delair",
        "fullExtent": get_mercator_full_extent(ortho_dataset)
    }

def create_arcgis_basemap(title, layers):
    return {
        "baseMap": {
            "baseMapLayers": [
                {
                    "id": "defaultBasemap",
                    "layerType": "ArcGISTiledMapServiceLayer",
                    "url": "https://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer",
                    "visibility": True,
                    "opacity": 1,
                    "title": "Topographic"
                }
            ],
            "title": title
        },
        "spatialReference": {
            "wkid": 102100,
            "latestWkid": 3857
        },
        "authoringApp": "WebMapViewer",
        "authoringAppVersion": "7.1",
        "version": "2.14",
        "operationalLayers": layers,
    }

def publish_map_on_arcgis(title, layers):
    arcgis_layers = [create_arcgis_layer(name=layer.get('name'), url=layer.get('url')) for layer in layers]
    base_map = create_arcgis_basemap(title, arcgis_layers)
    item_prop = {
        "type": "Web Map",
        "title": title,
        "tags": ["WebTiledLayer", "delair.ai"],
        "snippet": "",
        "text": base_map,
        "extent": [[bbox[0], bbox[1]],[bbox[2], bbox[3]]]
    }
    
    generated_map = gis.content.add(item_properties=item_prop)
    return generated_map

Let's generate a map on Arcgis with both Orthomosaic and DSM¶

In [ ]:

custom_dsm_tile_url = generate_custom_tile_url(
    dsm_tile_url,
    min=188,
    max=201,
    outofrange_type='transparent'
)

layers = [
    {'name': 'Orthomosaic', 'url': ortho_tile_url},
    {'name': 'DSM', 'url': custom_dsm_tile_url}
]

generated_map = publish_map_on_arcgis('Serving delair.ai data on Arcgis', layers)
generated_map

In [ ]:

m = gis.map(generated_map)
m

In [ ]:

Share delair.ai map tiles with the Python SDK¶

You may execute this notebook one cell after the other with Shift + Enter¶

Requirements¶

Download sample files (images, mesh, raster...)¶

Create the project¶

Create and upload the orthomosaic¶

DSM (Digital Surface Model)¶

After the raster uploads, they need to be tiled. It can take up to 5 minutes.¶

Tile layer URLs Generation¶

Determine the center of the raster and its bounding box¶

Visualization with Folium¶

Orthomosaic¶

With Argis¶

New map¶

Let's generate a map on Arcgis with both Orthomosaic and DSM¶

You may execute this notebook one cell after the other with `Shift + Enter`¶