Použití topografických dat SRTM (Česká přehrada)

4.4.2021 00:51 | Přečteno: 4498× | Linux | | poslední úprava: 4.4.2021 00:51

SRTM jsou výškopisná data (topografie) zemského povrchu - jsou tam vidět kopce a údolí. Ukážeme si, jak s nimi pracovat. A jako vedlejší efekt si ukážeme, jak vygenerovat pro webový prohlížeč vlastně jakoukoli mapu.

SRTM data se dají stáhnout například zde. Jejich rozlišení je 1 úhlová vteřina, tj. v našich zeměpisných šířkách obvod_Země/360/3600 = 31 m na výšku a *cos(50°) = 20 metrů na šířku. Soubory se jmenují například N50E014.SRTMGL1.hgt.zip a tento soubor obsahuje data od 50° do 51° severní šířky a 14° do 15° východní délky, což jsou cca. severní Čechy. Po stažení jej standardním způsobem rozzipujeme.

Soubor obsahuje matici 3601x3601 16bit intů a je tedy dlouhý 3601*3601*2 = 25934402 bajtů (26 MB). Poněkud překvapivé je, že čísla jsou big endian. Čísla přímo vyjadřují nadmořskou výšku v metrech (většina rozsahu je tedy nevyužita, neboť nic ani zdaleka se blížící 65535 metrům se na Zemi nevyskytuje). Soubor můžeme načíst v Pythonu a zobrazit:

import numpy as np
import matplotlib.pyplot as plt

data = np.fromfile("N50E014.hgt", dtype=np.uint16)
data = data.reshape((3601, 3601))
data.byteswap(inplace=True)
plt.imshow(data, cmap="turbo")
plt.show()

Vidíme tam třeba údolí různých řek, hurá. Vidíme také (pokud víte jak mají severní Čechy vypadat), že obrázek je „zplácnutý“, protože pixely nejsou čtvercové, ale 30x20 metrů. Obrázek vyexportujeme, a tentokrát použijeme barevnou škálu terrain, která se k tomuto hodí:

def export_figure_matplotlib(arr, f_name, dpi=200):
    fig = plt.figure(frameon=False)
    fig.set_size_inches(arr.shape[1]/dpi, arr.shape[0]/dpi)
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)
    ax.imshow(arr, cmap="terrain")
    plt.savefig(f_name, dpi=dpi)

export_figure_matplotlib(data, "a.png")

(ano, je to příšerné, a jak exportuju obrázky normálně popíšu v příštím článku)

Nyní obrázek plácneme do leafletu. Jednoduše nastavíme levý horní a pravý dolní okraj a ono ho to tam nějak prskne.

<html>
<head>
 <link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css"
   integrity="sha512-xodZBNTC5n17Xt2atTPuE1HxjVMSvLVW9ocqUKLsCC5CXdbqCmblAshOMAS6/keqq/sMZMZ19scR4PsZChSR7A=="
   crossorigin=""/>
 <script src="https://unpkg.com/leaflet@1.7.1/dist/leaflet.js"
   integrity="sha512-XQoYMqMTK8LvdxXYG3nZ448hOEQiglfqkJs1NOQV44cWnUrBc8PkAOcXy20w0vlaXaVUearIOBhiXZ5V3ynxwA=="
   crossorigin=""></script>
</head>
<body>
<div id="mapid" style="width: 900px; height: 800px;"></div>
<script>
var map = L.map('mapid').setView([50.5, 14.5], 8);
L.tileLayer('https://mapserver.mapy.cz/turist-m/{z}-{x}-{y}', {
  attribution: 'Map data &copy; <a href="https://www.seznam.cz">Seznam.cz, a.s.</a>'
}).addTo(map);
img = L.imageOverlay("a.png", [[51, 14], [50,15]], {opacity: 0.8, zIndex:20}).addTo(map);
</script>
</body>
</html>

Jako zdroj dlaždic jsem použil mapy.cz, protože jsou prostě nejlepší (druhé nejlepší jsou oficiální dlaždice openstreetmap, ale ty často ratelimitují k nepoužitelnosti). Blbé je, že licence mluví jenom o použití s jejich vlastním javascriptovým zobrazovačem, a s Leafletem se to nejspíš nesmí (na druhou stranu copyright uvedený je a ničím jiným se to neliší). Osobně jsem si alespoň spustil lokální kešovací proxy ať jim při svém používání nevytěžuju server. Zde jsem jako podklad použil turist-m (turistická mapa), ale pro některé aplikace je zajímavá hybrid-base-m, což je poloprůhledná mapa s kontrastními významnými liniemi (hranice, města, silnice…), přes kterou se dá dobře kreslit vlastní obsah (je to to, co se na mapy.cz ukazuje jako překryv letecké pokud si nedáte čistou leteckou z menu).

Vidíme, že to docela sedí na featury v mapě. Nevím, jestli to má sedět úplně, nebo je tohle špatná projekce (vlastně jsme jenom natáhli obdélník aby pasoval v rozích), každopádně tohle natažení obdélníku používám i odvážněji (data z radaru nakreslím na placku a pak ji takhle narvu do mapy, a to je zaručeně špatně, protože zakřivení Země) a do vzdálenosti cca. 150 km (tj. obdélník 300x300 km) to funguje bez nějakého zjevného nepasování.

Zpracování větší plochy

Pokud chceme stavět rozsáhlejší vodní díla, budeme potřebovat spojit více těchto dlaždic 1°x1° do jednoho obrázku. Samozřejmě můžeme nějak lepit ta pole, ale zajímavé bude zkusit použít gdal, čímž získáme rovnou obrázky georeferencované a vůbec. Vytvoříme virtuální dataset, což je jakoby veliký obrázek, ale ve skutečnosti malé XML, které popisuje, že pro souřadnice X a Y se má sáhnout do nějakého souboru na disku, a pro jiné souřadnice do jiného souboru.

Druhá věc, kterou zde používám, je falešná cesta začínající /vsigzip/ - já jsem si soubory po vybalení ze zipu zagzipoval, aby nebyly moc velké, a tohle je „driver“, který umožňuje gdalu transparentně přistupovat k zagzipovaným souborům.

$ gdalbuildvrt /home/jenda/brmdam/all.vrt /vsigzip//home/jenda/simulator/data/N43E015.hgt.gz /vsigzip//home/jenda/simulator/data/N43E016.hgt.gz [a další]
$ cat /home/jenda/brmdam/all.vrt
<VRTDataset rasterXSize="36001" rasterYSize="25201">
...
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/vsigzip//home/jenda/simulator/data/N43E015.hgt.gz</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="3601" RasterYSize="3601" DataType="Int16" BlockXSize="3601" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="3601" ySize="3601" />
      <DstRect xOff="3600" yOff="21600" xSize="3601" ySize="3601" />
      <NODATA>-32768</NODATA>
    </ComplexSource>

Takhle vytvořený dataset můžeme třeba převést do geotiffu.

$ gdalwarp
    -r bilinear # použitá interpolace
    -t_srs EPSG:3857 # cílová projekce - toto používají online mapy
    -ts 2048 0 # cílová velikost v pixelech (0 = druhou souřadnici dopočítá)
    ~/brmdam/all.vrt /tmp/all.tiff

Takto vytvořený TIFF můžeme otevřít třeba v GIMPu (doporučuji Colors → Levels, protože jinak je to šedé a nic tam neuvidíte), ale zajímavé je, že to není jenom obyčejná bitmapa, ale že to je georeferencované:

$ gdalinfo /tmp/all.tiff
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
...
Upper Left  ( 1558457.410, 6446299.894) ( 13d59'59.50"E, 50d 0' 0.50"N)
...

Díky tomu můžeme používat různé další tooly co gdal a osgeo nabízí, třeba z toho udělat dlaždice pro zobrazení online mapy v různém zoomu:

$ gdal2tiles.py /tmp/all.tiff tile/

Ty pak můžeme přímo načíst jako vrstvu do Leafletu:

var tiles = L.tileLayer('./tile/{z}/{x}/{y}.png', {tms: true, opacity: 0.7, minZoom: 5, maxZoom: 8}).addTo(map);

Tohle je super, protože návštěvníci vašeho webu nebudou tahat gigapixelový obrázek, ale budou mít normální dlaždicovou mapu. A navíc to má korektně projekci.

Kreslení vlastních dat do mapy

Naším cílem ale je umět si do mapy s výškopisem kreslit vlastní věci. To uděláme tak, že místo tiffu uděláme PPM, a v něm pak vyměníme obrazová data (přičemž to celé zůstane georeferencované). Využijeme zde triku, že když se ukládá do formátu, který nemá standardizovaná geo metadata, tak je gdal uloží bokem do název_souboru.aux.xml. No a ten druhý soubor je obyčejný obrázek, se kterým můžeme manipulovat.

$ gdalwarp -r bilinear -t_srs EPSG:3857 
  -ot Byte # PPM neumí uložit 16bit data, takže je musíme nějak naškálovat, ale to je jedno, protože si obrázek stejně nahradíme vlastním
  -ts 2048 0 ~/brmdam/all.vrt /tmp/all.ppm

# vznikl též soubor /tmp/all.ppm.aux.xml

Takže teď máme GeoTIFF a PPM a pixely na obou jsou na stejných „pozemských“ souřadnicích. GeoTIFF jde normálně otevřít PILem a načtená čísla stále odpovídají nadmořské výšce v metrech:

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

im = Image.open("/tmp/all.tiff")
data = np.array(im)

# nedefinovaná data (moře) jsou -2^15 což rozbije autoscaling barviček
# proto nastavíme maximální nadmořskou výšku na 0
data[data<0] = 0

plt.imshow(data)
plt.show()

Z hodnot můžeme třeba ručně udělat RGB hodnoty color mapy terrain a výsledek uložit, tentokrát ale do PNG, protože gdal2tiles.py neumí načíst PPM (gdalwarp zase pro změnu neumí zapsat PNG).

from matplotlib import cm

# funkce cm.colormap žere buď inty od 0 do 255 nebo floaty od 0 do 1 a vrací RGBA barvu
# pro intový vstup mapuje [0..255] → [0..255]^4
# pro floatový vstup mapuje [0..1] → [0..1]^4
# v tomhle se snadno udělá chyba…
# my budeme používat floaty, mapu normalizujeme od 0 do 1

data_normalized = data / np.amax(data)

data_colorized = cm.terrain(data_normalized)
plt.imshow(data_colorized)
plt.show()

# nakreslíme někam doprostřed červený obdélník protože můžeme
data_colorized[1000:1100,1000:1200,:] = [1,0,0,1]

im = Image.fromarray((data_colorized*255).astype(np.uint8))
im.save("/tmp/all.png")

$ cp /tmp/all.ppm.aux.xml /tmp/all.png.aux.xml

Tím jsme získali georeferencované PNG s naším vlastním obsahem a z něj můžeme vyrobit třeba zase ty dlaždice (všechny nástroje spolu s obrázkem automaticky načtou i příslušný .aux.xml soubor).

$ gdal2tiles.py /tmp/all.png tile2/

Zaplavování

Potřebujeme umět vyplnit oblast stavěné přehrady. Myslel jsem si, že to je jednoduché omezené DFS na pár řádků, ale nepovedlo se mi to napsat tak, aby mi u velkých obrázků nedošel stack a vůbec je to celé nějaké složité a FFFUUU. Naštěstí existuje floodfill v skimage.

Nejdřív si ale musíme nadefinovat hráze. To jsem udělal tak, že jsem v GIMPu nakreslil bílou hráz a všude okolo černo. Obrázek pak můžeme načíst a hodnoty použít.

# to PNG je RGBA (ano, jde to nastavit při exportu z GIMPu) tak z toho vezmeme jenom první složku
hraze = np.array(Image.open("/tmp/hraze.png"))[:,:,0]

# budeme napouštět do nadmořské výšky 200 m
hladina = 200

# nastavíme hráze nad úroveň hladiny
data[hraze>0] = hladina+1

# nastavíme co potenciálně může být pod vodou na -1
data_flooded = data.copy()
data_flooded[data<hladina] = -1

# spustíme zaplavování, které nastaví zaplavené oblasti na -2
# jako počáteční bod použijeme ten obdélník z minula
from skimage.segmentation import flood_fill
result = flood_fill(data_flooded, (1000, 1000), -2)

# promítneme zaplavené oblasti do původního obrázku
# nastavíme tam -400 aby byla zřetelnější hranice mezi přehradou a břehem
data[result == -2] = -400
data += 400

# normalizujeme a obarvíme
data_normalized = data / np.amax(data)
data_colorized = cm.terrain(data_normalized)

# hráze nakreslíme červeně
data_colorized[hraze>0] = [1,0,0,1]

# uložíme
im = Image.fromarray((data_colorized*255).astype(np.uint8))
im.save("/tmp/all.png")

Poslední věc, kterou vidíte na ceskaprehrada.cz, je speciální barevná mapa pro zaplavené (různé odstíny modré znázorňující hloubku) a nezaplavené (zelená) oblasti, kterou se mi už nechce vysvětlovat, ale předpokládám, že je to zjevné.

Komprese dat

HGT soubory lze zmenšit na devítinu (výsledek bude mít 1201x1201 pixelů), gdal to automaticky detekuje podle velikosti souboru a vše bude fungovat. To se hodí pokud je zpracování větší oblasti příliš pomalé. Tady je program (součást většího bazmeku který chceme někdy uvolnit), který to dělá s maximovým filtrem.

#include <stdlib.h>
#include <stdio.h>
#include <inttypes.h>
#include <stdbool.h>
#include <unistd.h>
#include <arpa/inet.h>

int16_t* inbuf;

#define INBYTES (3601*3601*2)

#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))
#define CLAMP(x, lower, upper) (MIN(upper, MAX(x, lower)))
#define PE(...) do{ fprintf( stderr, __VA_ARGS__ ); } while( false )

int main(int argc, char** argv) {
//    PE("Downsample a SRTM HGT file with 30-meters resolution using 3x3 maxpool\n");
//    PE("Usage: %s\n", argv[0]);

  inbuf = malloc(INBYTES);
  if(!inbuf) {
    perror("malloc");
    return 1;
  }
  size_t r = fread(inbuf, 1, INBYTES, stdin);
  if(r != INBYTES) {
    perror("read");
    return 1;
  }
  for(int j = 0; j<1201; j++) {
    for(int i = 0; i<1201; i++) {
      int sx = i*3-(i/401);
      int sy = j*3-(j/401);
      //if(sx == 3) printf("%i -> %i\n", j, sy);
      uint16_t res = 0;
      for(int ii = sx; ii<sx+3; ii++) {
        for(int jj = sy; jj<sy+3; jj++) {
          //printf("%i %i %i %i\n", i, j, ii, jj);
          res = MAX(res, ntohs(inbuf[jj*3601+ii]));
        }
      }
      res = htons(res);
      fwrite(&res, 2, 1, stdout);
    }
  }
  free(inbuf);
}

Další služby využívající SRTM

BRouter je plánovač cest, který zohledňuje kopce.

Odkazy

• Pannonian Sea
• Kobylské jezero
• Čejčské jezero

       

Komentáře

[_]

Od čeho se odráží ten radar? Od povrchu země nebo od vegetace? Nebo je přesnost tak malá, že to nemá smysl řešit (i když stromy mohou být vysoké)?

[Jendа]

Tak nějak od všeho. A stromy stíní, a na několika našich lokacích je to problém.

[Radovan]

To mi připomnělo: http://technomorous.eu/post/173601337433 Plus předchozí a následující post.

[Jendа]

To je divné že tam píše že to zpracovával tak dlouho, všechny skripty v zápisku doběhnou za řádově sekundy, maximálně desítky sekund, a to není nijak optimalizované.

[Hx.]

Ne každý má tak rychlý počítač jako ty, Jirko.

[plostenka]

Ten podivny rozmer 1021x1201 (NASA SRTM) nebo 3601x3601 (ESA SRTM) je zpusobeny tim, ze posledni vychodni sloupec, respektive jizni radek je zduplikovany se sousedni dlazdici. Vyskova nepresnost je cca 10 metru a chyba mereni je zanesena jako 2^16-1 (civilni SRTM je unsigned nad hladinou more).

[❗Červivý Hřib krade❗]

zase jako tu velikonoční voblejvačku neber tak vážně :D :D ;D ;D

btw proč je vtomdletom

#define PE(...) do{ fprintf( stderr, __VA_ARGS__ ); } while( false )

jakoby to do while false?????? :O :O

[Jendа]

Je to idiom jak pomocí define vyrobit něco, co se bude chovat jako funkce (tj. půjde napsat PE(blabla);).

[❗Červivý Hřib krade❗]

supr dík ;D

[Dirka]

Mozna by se ti mohla libit https://www.freemap.sk/?map=14/50.851366/15.227480&layers=X (layer X mi pride nejlepsi pro vrstevnice, je OSM based).

Pouzivani dat z osm.org se nedoporucuje, vhodne je rozbehnout si vlastni dlazdice.

Ty oblasti by sly asi i vycitat pomoci OSM dat, tzn u vodniho dila to bude relace.

Klidne se ptej, pokud te neco ohledne OSM zajima.

[Člověk z Horní Dolní]

Existuje nebo se chystá něco lepšího než SRTM, ale stále zdarma?

[Martin Tůma]

Pro mnohé evropské země jsou k dispozici LiDAR data. Stejně jako u map se už dneska dost často dají volně použít data poskytovaná přímo národníma "geo-agenturama".

[...]

to je taky pekne hruba sit

[Martin Tůma]

Rozlišení to má stejný (1") jako SRTM, ale obzvláště v Alpskejch regionech jsou ty data mnohem přesnější. Na běžný "domácí hraní" si s GPS logama je to zcela dostačující, ale jestli podle toho chceš stavět aquadukt, tak je to samozřejmě málo...

[...]

chci na tom vidět/změřit výšky hor a sedel, ale to jsou struktury výrazně menší

[Senior Fullstack Javascript Developer]

Na Slovensku robí letecké laserové zameriavanie Úrad geodézie, kartografie a katastra SR (ÚGKK SR) určite máte nejaký taký projekt aj v ČR. Zatiaľ majú polovicu Slovenska a budú to robiť ěste asi dva roky:

https://www.geoportal.sk/sk/udaje/lls-dmr/?nastavene=1

Vyskova presnosť je v milimetroch: https://www.geoportal.sk/files/zbgis/lls/parametre-lokality-zberu-udajov-lls.pdf Napríklad vo vojenskom cvičnom priestore vidieť diery po rôznej munícii: https://www.freemap.sk/?map=17/48.479180/17.123464&layers=X Alebo rôzna drobná historická povrchová ťažba (tzv. pingy): https://www.freemap.sk/?map=17/48.783713/18.069892&layers=XI

[j]

Rek bych, ze mnohem realnejsi scenar bude, kdyz nechas ujet ten ujizdejici svah pod D8 = hraz vybudujes v udoli labe.

Prirozenej prepad takovy prehrady je nekde v okoli Mostu, zpet by to pak teklo korytem Biliny. A Praha by konecne mela svoje more, protoze nekde na okraji by byl breh (teda samo, centralni cast by byla spolehlive pod vodou).

Bywoko - hladina vltavy v praglu cca 180m, Labe v prackovicich cca 140m, zakladni vyska toho sesuvu rekneme 400m. Rekneme ze by to zvladlo hraz o vysce 100m => naprosto spolehlive by to nepreteklo tady.

To misto kudy by to teklo ma neco kolem 200m. Takze by v praglu v udoli bylo cca 20m vody ;D

---

Dete s tim guuglem dopice!

[Ondřej J]

Mě by zajímalo, jak dlouho by se to napouštělo (Arvíz). Asi při zanedbání výparu, protože s výparem by se to imho na cílových 200m nenapustilo nikdy.