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

[[Category:Croatia/Kratki_vodici|D]]

= Digitalni model terena iz CSV datoteke pomoću GDALa =

Kako bi uspješno pratili kratki vodič morate zadovoljiti sljedeće:
# instalirati GDAL programski paket za vaš operativni sustav
# preuzeti i otpakirati [http://hr.osgeo.org/sites/default/files/podaci.zip ZIP datoteku] koja sadrži podatke u CSV formatu

Datoteka ''podaci.csv'' sadrži sljedeće zapise:
<pre>
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
...
</pre>

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

Tijek procesa:
# napraviti raster pomoću ''gdal_grid''
## definirati CSV kao virtualni format za OGR
## odabranom metodom interpolacije napraviti raster
# izraditi slojnice pomoću ''gdal_contour''

 
== ''Gridding'' pomoću gdal_grid ==

=== VRT virtualni format ===

[http://www.gdal.org/gdal_vrttut.html VRT] je tzv. virtualni format, odnosno tekstualna datoteka s metapodacima koji opisuju neki set podataka.

[[Image:Points.png|thumbnail|frameless| Točke sa visinskim vrijednostima]] 

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

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:
<pre>
$ 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)
</pre>

Primjerice znamo da skup podataka ima 4448 točaka i da nema definiran prostorni referentni sustav (SRS). SRS možemo definirati tako da dodamo:
<pre>
<LayerSRS>+proj=tmerc +ellps=bessel +k_o=0.9999 +x_0=5500000 +lon_0=15 +units=m</LayerSRS>
</pre>
č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 [http://www.gdal.org/gdal_grid.html 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: [[Image:Dem.png|thumbnail|frameless| Interpolacijom dobiveni raster DEM]] [[Image:Interpolation.png|thumbnail|frameless| Interpolacijom dobiveni raster DEM s prikazom tocaka]] 
;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:
<pre>
gdal_grid -a invdist:power=2.0:smoothing=1.0:radius1=20:radius2=20 -of GTiff -ot Float32 -l podaci podaci.vrt dem.tiff
</pre>

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:
<pre>
$ 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
</pre>

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 <tt>-outsize 319 466</tt> prije navedenoj ''gdal_grid'' naredbi.

== ''gdal_contour'' izrada slojnica ==
[[Image:Slojnice.png|thumbnail|frameless| Slojnice dobivene sa gdal_contour alatom]] 
Alat za izradu slojnica [http://www.gdal.org/gdal_contour.html gdal_contour] omogućava spremanje slojnica u bilo kojem vektorskom formatu koji podržava OGR, no ukoliko ne definiramo format spremiti će ih u ''SHP'' formatu.

Naredbom:  
<pre>
gdal_contour -a visina -3d -i 1 dem.tiff slojnice1.shp
</pre>

se u tekućem direktoriju stvara datoteka ''slojnice1.shp'' koristeći ''dem.tiff'' izvor podataka sa intervalom od '''1''' jedinice (u ovom slučaju metra). Osim toga svaka točka će biti određena u 3d koordinatnom sustavu, a informacija o visini će se zapisati kao atribut ''-a visina''. 

Vizualizacija ovisi o softveru koji preferirate. Slike u tutorijalu su napravljene vizualizacijon deriviranih podataka u QGIS alatu.
[[Image:iso.png|center|thumbnail|frameless| Vizualizacija svih slojeva dobivenih geoprocesiranjem opisanim u tutorijalu]]