Using `geom_imshow()` to draw DEM on map¶

SRTM1N00E012V3 courtesy of the U.S. Geological Survey

In [1]:

import numpy as np
from PIL import Image

from lets_plot import *
from lets_plot.geo_data import *
LetsPlot.setup_html()

The geodata is provided by © OpenStreetMap contributors and is made available here under the Open Database License (ODbL).

In [2]:

def get_dem_data(url):
    import requests
    from osgeo import gdal
    r = requests.get(url)
    vsipath = '/vsimem/dem'
    gdal.FileFromMemBuffer(vsipath, r.content)
    return gdal.Open(vsipath)

In [3]:

plot_size = 600
subplot_size = 500

In [4]:

# Data from geocoding
states = geocode_states().scope("Gabon").inc_res(4).get_boundaries()
cities = geocode_cities().scope("Gabon").get_centroids()

In [5]:

# Data from DEM
dataset = get_dem_data("https://raw.githubusercontent.com/JetBrains/lets-plot-docs/master/data/ivindo_river.tif")
dataset_array = dataset.ReadAsArray()
transform = dataset.GetGeoTransform()
lon_min = transform[0]
lon_max = lon_min + (transform[1] * dataset_array.shape[1])
lat_max = transform[3]
lat_min = lat_max + (transform[5] * dataset_array.shape[0])

In [6]:

plot = ggplot() + \
    geom_imshow(dataset_array, extent=[lon_min, lon_max, lat_min, lat_max], cmap="viridis") + \
    geom_map(data=states, color="black", alpha=0) + \
    geom_point(data=cities, shape=21, size=4, color="black", fill="white", \
               tooltips=layer_tooltips().line("@{found name}")) + \
    ggsize(plot_size, plot_size) + \
    theme_classic() + theme(axis="blank")
plot + ggtitle("Ivindo river in Gabon")

Out[6]:

Using geom_imshow() to draw DEM on map¶

Using `geom_imshow()` to draw DEM on map¶