Croatia/Kratki vodici/digitalni-model-terena-iz-csv-datoteke-pomocu-gdala

Digitalni model terena iz CSV datoteke pomoću GDALa

Kako bi uspješno pratili kratki vodič morate zadovoljiti sljedeće:

1. instalirati GDAL programski paket za vaš operativni sustav
2. preuzeti i otpakirati ZIP datoteku koja sadrži podatke u CSV formatu

Datoteka podaci.csv sadrži sljedeće zapise:

X,Y,Z
5585884.50,5063582.05,63.45
5585884.67,5063582.95,63.90
5585884.75,5063583.81,63.97
5585884.65,5063584.70,64.18
5585884.27,5063585.74,64.30
5585883.73,5063586.89,64.37
...

Zaglavlje datoteke definira koje informacije sadrži koja kolona, u ovom slučaju X i Y koordinatu te visinu Z.

Tijek procesa:

1. napraviti raster pomoću gdal_grid
   1. definirati CSV kao virtualni format za OGR
   2. odabranom metodom interpolacije napraviti raster
2. izraditi slojnice pomoću gdal_contour

Gridding pomoću gdal_grid

VRT virtualni format

VRT je tzv. virtualni format, odnosno tekstualna datoteka s metapodacima koji opisuju neki set podataka.

<OGRVRTDataSource>
  <OGRVRTLayer name="podaci">
     <SrcDataSource>podaci.csv</SrcDataSource> 
     <GeometryType>wkbPoint</GeometryType>
     <GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/> 
  </OGRVRTLayer>
</OGRVRTDataSource>

Zapis <OGRVRTLayer name="podaci"> definira naziv sloja podataka, <GeometryType> tip podataka, koji je u ovom slučaju točka, i <GeometryField koji je definiran zapis, odnosno u kojoj koloni se nalazi koji podatak. Naziv sloja podataka bi trebao biti isti kao naziv datoteke bez ekstenzije, dakle podaci i podaci.

Tekstualnu datoteku možemo nazvati kako želimo, ali za potrebe ovog vodiča neka bude podaci.vrt.

Nakon spremanja datoteke možemo iskoristiti alat ogrinfo kako bi saznali osnovne informacije o ovo skupu podataka:

$ ogrinfo podaci.vrt podaci -so
INFO: Open of `podaci.vrt'
      using driver `VRT' successful.

Layer name: podaci
Geometry: Point
Feature Count: 4448
Layer SRS WKT:
(unknown)
X: String (0.0)
Y: String (0.0)
Z: String (0.0)

Primjerice znamo da skup podataka ima 4448 točaka i da nema definiran prostorni referentni sustav (SRS). SRS možemo definirati tako da dodamo:

<LayerSRS>+proj=tmerc +ellps=bessel +k_o=0.9999 +x_0=5500000 +lon_0=15 +units=m</LayerSRS>

čime bi definirali 5 zonu GaussKruger projekcije bez transformacijskih parametara prema wgs84 u PROJ.4 formatu.

Gridding i odabir metode interploacije

Već je spomenuto kako gdal_grid podržava tri metode interpolacije (link). Osim toga za svaku se metodu mogu kontrolirati različiti parametri. Najkorisnije je spomenuti da bez definicije parametara gdal_grid u obzir uzima sve točke zadanog područja, što vrlo vjerojatno nije željno ponašanje.

Od parametara su najzanimljiviji:

radius1

prvi promjer elipse područja traženja, ukoliko je 0 u obzir se uzima cijelo područje

radius2

drugi promjer elipse područja traženja, ukoliko je 0 u obzir se uzima cijelo područje

angle

kut zakreta elipse u smjeru kazaljke na satu

min_points

minimalan potreban broj točaka u traženom području, ukoliko ima manje točaka od traženog broja, ta točka se smatra bez vrijednosti nodata

nodata

vrijednost koja predstavlja 'nepostojeću vrijednost'

Recimo da želimo koristiti metodu invdist potpuna naredba bi izgledala:

gdal_grid -a invdist:power=2.0:smoothing=1.0:radius1=20:radius2=20 -of GTiff -ot Float32 -l podaci podaci.vrt dem.tiff

Ukratko, postavili smo parametre metode interpolacije, kao izlazni format -of GTiff postavljen je GeoTiff, tip izlazne datoteke je -ot Float32 što znači da je svaki element rastera 32bitni realni broj (i više nego dovoljno za naše potrebe). Također smo odabrali sloj podataka -l podaci, postavili ulaznu datoteku podaci.vrt i izlaznu datoteku dem.tiff.

Nakon izvršavanja naredbe upotrebom alata gdalinfo možemo saznati osnovne informacije o novo stvorenom rasteru:

$ gdalinfo dem.tiff
Driver: GTiff/GeoTIFF
Files: dem.tiff
Size is 256, 256
Coordinate System is `'
Origin = (5585728.379999999888241,5063425.360000000335276)
Pixel Size = (1.239687500001310,1.811992187500437)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 5585728.380, 5063425.360) 
Lower Left  ( 5585728.380, 5063889.230) 
Upper Right ( 5586045.740, 5063425.360) 
Lower Right ( 5586045.740, 5063889.230) 
Center      ( 5585887.060, 5063657.295) 
Band 1 Block=256x8 Type=Float32, ColorInterp=Gray

Između ostalog i rezoluciju pojedinog elementa rastera, ukoliko želimo raster drugačije rezolucije potrebno je za gdal_grid postaviti parametar -outsize. Primjerice za rezoluciju 1x1 dodali bi parametar -outsize 319 466 prije navedenoj gdal_grid naredbi.