Difference between revisions of "Croatia/Kratki vodici/digitalni-model-terena-iz-csv-datoteke-pomocu-gdala"
Wiki-Dodobas (talk | contribs) |
Wiki-Darko (talk | contribs) |
||
(6 intermediate revisions by one other user not shown) | |||
Line 1: | Line 1: | ||
+ | [[Category:Croatia/Kratki_vodici|D]] | ||
+ | |||
= Digitalni model terena iz CSV datoteke pomoću GDALa = | = 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]] |
Latest revision as of 11:31, 11 December 2010
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 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:
- 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
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.
gdal_contour izrada slojnica
Alat za izradu slojnica 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:
gdal_contour -a visina -3d -i 1 dem.tiff slojnice1.shp
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.