Manipularea datelor spațiale folosind GDAL. I Date raster

Introducere

GDAL (GDAL – Geospatial Data Abstraction Library) este o librarie opensource pentru manipularea datelor geospatiale, atat in format raster cat si vector. Aceasta este compusa dintr-o suita de tooluri command line, pe care le vom prezenta in acest articol, si din cateva API-uri conexe ce permit dezvoltarea unor scripturi/tooluri in orice limbaj de comoditate Perl, Python, Java, sau limbaje precum C / C++ / C#.

Libraria suporta orice platforma Uni(x), precum si platforme Win32 sau Mac OS. GDAL este publicat sub licenta Public Domain, codul sursa sau distributiile binare pentru Linux si Win32 putand fi descarcat de pe site-ul http://www.gdal.org.

Suita GDAL este considerata la ora actuala cea mai puternica si vasta suita de manipulare raster/vector, stand la baza multor softuri proprietare, fiind adoptat si integrat in solutiile GIS ale unor companii foarte cunoscute precum: ESRI, Leica, GlobalMapper, sau in celebrul Google Earth.

GDAL se detaseaza net de orice adversar prin numarul mare de formate raster si vector suportate. Cele mai importante formate de date spatiale, precum GeoTiff, ERDAS Img, ESRI Shapefile, sunt foarte bine implementate si optimizate la manipulare. GDAL implementeaza suport solid si extensiv pentru formate non-spatiale, cum sunt formatele si derivatele posibile de rasteri gen BMP, PNG, JPEG, TIFF. Totodata, GDAL implementeaza sintagma SQL atat pentru metadatele (ex: .dbf, .tab) ce vin atasate unor formate vectoriale cat si suport pentru conexiuni directe la servere SQL, atat opensource cat si proprietare. GDAL mai suporta manipularea directa a materialelor disponibile online via map servere ce implementeaza protocolul WMS.

Modulului de transformare a coordonatelor utilizeaza baza de date EPSG (http://www.espg.org), disponibila si aceasta sub o licenta gratuia, cu anumite restrictii, perfect de inteles, la capitolul re-distribuire/alterare. Baza EPSG cuprinde un numar impresionant de sisteme de coordonate, inclusiv cele romanesti, orice profesionist putand regasi aici toate transformarile si datum-urile necesare.

O lista detailata a formatelor spatiale raster suportate de GDAL se gaseste la http://www.gdal.org/formats_list.html. Lista similara pentru cele vectoriale poate fi accesata la http://www.gdal.org/ogr/ogr_formats.html. Site-ul contine si documentatia integrala GDAL (mai putin exemplele de utilizare).

Suita GDAL este inclusa in distributii linux precum Fedora, Debian sau Gentoo, impreuna cu librariile EPSG. Mai mult, suita este introdusa in path-urile de executie pentru a usura munca celuia care il folseste si a-l scuti de instalari anevoioase si extensive. Manualele aferente sunt disponibile prin celebra comanda „man”, nefiind necesara accesarea on-line site-urilor GDAL sau EPSG.

Pentru platforme Win32 exista un re-pack, intitulat FWTools si intocmit de insusi autorul GDAL: Frank Warmerdam, disponibil gratuit pentru download pe site-ul http://fwtools.maptools.com. FWTools functioneaza impecabil in conjunctie cu mediul Cygwin (http://cygwin.com), o consola foarte utila pentru cei care nu vor sa renunte la platforma Windows si doresc in acelasi timp si un mediu de scripting pentru interpretorul de comenzi „bash”, intr-o fereastra scalabila la orice dimensiuni, net superior consolei DOS disponibila nativ in Win32.

In acest prim articol o sa incercam sa introducem doar partea de raster, urmand ca in articolele viitoare sa prezentam extensiv si partea vectoriala, respectiv partea de integrare in scripturi de uz general.

GDAL in actiune

Vom incepe prin a descrie exact ce contine suita GDAL. Din punct de vedere raster avem:

• gdalinfo – afiseaza informatii despre un obiect/fisier raster;
• gdal_translate – converteste obiecte raster dintr-un format de fisier in altul;
• gdaladdo – adauga preview/piramida de previzualizare la un raster;
• gdalwarp – converteste un obiect raster dintr-o proiectie in alta;
• gdaltindex – genereaza un index al subseturilor pentru utilizarea in Mapserver;
• gdal_contour – extrage curbe nivel din obiecte raster de elevatie;
• rgb2pct.py – converteste un obiect raster din 24 biti RGB in 8 biti indexat;
• pct2rgb.py – converteste un obiect raster din 8 biti indexat in 24 biti RGB;
• gdal_merge.py – combina mai multe obiecte raster intr-un singur obiect raster;
• gdal2tiles.py – creeaza o structura TMS pentru vizualizare web KML (tip google);
• gdal_rasterize – rasterizeaza vectori intr-un singur obiect/fisier raster;
• gdaltransform – transforma coordonate dintr-un sistem in altul;
• nearblack – converteste pixelii frontiera cu nuante apropiate de alb sau negru in valori exacte alb/negru;
• gdal_retile.py – reconstruieste piramidele pentru un set de subseturi (tile-uri);
• gdal_grid – interpoleaza seturi de date dispersate pentru a obtine un grid;
• gdal-config – obtine diverse informatii referitoare la libraria GDAL instalata (disponibila doar pe platforme Unix);

Din punct de vedere vectorial avem:

• ogrinfo – afiseaza detalii despre obiectele vectoriale;
• ogr2ogr – converteste/transforma/reproiecteaza obiecte vectoriale;
• ogrtindex – creeaza poligoane vectoriale dupa extendul spatial al obiectelor raster.

In continuare vom folosi notiunea de obiect pentru a face referire catre un set de date raster sau vector. Consideram ca termenul obiect este mai corect decat cel de fisier, deoarece, de multe ori, un strat este format dintr-un set de fisiere (ex: tif+tfw sau shp+shx+dbf). Mai mult, termenul fisier nu acopera satisfacator cazul staturile publicate via WMS (un strat al unui mapserver de examplu). In concluzie, notiunea de obiect se refera la un fisier sau grup de fisiere ce compun rasterul/vectorii, totoadata se refera la posibilitatea de fi local pe hardisk sau remote via http:// sau chiar WMS.

gdalinfo

Este un utilitar foarte prietenos ce ne permite sa aflam toate informatiile despre un fisier/obiect raster, idiferent daca acesta este fisier local sau este remote via http sau WMS. gdalinfo, ca si celelalte utilitare GDAL, face intodeauna abstractie de locatie si executa in mod similar input-ul.

In continuare va prezentam un exemplu de utilizare gdalinfo pe un obiect aflat remote:

$ gdalinfo /images/austrians/p_39_47_100.jpg
Driver: JPEG/JPEG JFIF
Files: none associated
Size is 3729, 5139
Coordinate System is `'
Image Structure Metadata:
INTERLEAVE=PIXEL
COMPRESSION=JPEG
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 5139.0)
Upper Right ( 3729.0,    0.0)
Lower Right ( 3729.0, 5139.0)
Center      ( 1864.5, 2569.5)
Band 1 Block=3729×1 Type=Byte, ColorInterp=Red
Image Structure Metadata:
COMPRESSION=JPEG
Band 2 Block=3729×1 Type=Byte, ColorInterp=Green
Image Structure Metadata:
COMPRESSION=JPEG
Band 3 Block=3729×1 Type=Byte, ColorInterp=Blue
Image Structure Metadata:
COMPRESSION=JPEG

Din rezultat aflam ca e vorba de un obiect raster, fara proiectie specificata, in trei benzi RGB, cu compresia de tip JPG si dimensiune de 3729×5139 pixeli. Pentru interpretare GDAL a folosit driverul generic raster JPEG (jfif). Practic acest obiect este o poza ne-georeferentiata.

Sa vedem un exemplu de obiect spatial adevarat, copiat local pe HDD de pe site-ul geo-spatial.org (mai exact linkul download/datele-landsat-etm-in-stereo701):

$ gdalinfo RGB321-L-34-020.tif 
Driver: GTiff/GeoTIFF
Files: RGB321-L-34-020.tif
Size is 2699, 2659
Coordinate System is:
LOCAL_CS[” Geocoding information not available Projection Name = Unknown Units = other GeoTIFF Units = other”,
UNIT[“unknown”,1]]
Origin = (235438.213641703130000,691205.184220034980000)
Pixel Size = (14.590382532428649,-14.590382532428643)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=IMAGINE TIFF Support
Copyright 1991 – 1999 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: etif.c $ $Revision: 1.10.1.9.1.9.2.11 $ $Date: 2004/09/15 18:42:01EDT $
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
TIFFTAG_RESOLUTIONUNIT=1 (unitless)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  235438.214,  691205.184) 
Lower Left  (  235438.214,  652409.357) 
Upper Right (  274817.656,  691205.184) 
Lower Right (  274817.656,  652409.357) 
Center      (  255127.935,  671807.271) 
Band 1 Block=64×64 Type=Byte, ColorInterp=Red
Band 2 Block=64×64 Type=Byte, ColorInterp=Green
Band 3 Block=64×64 Type=Byte, ColorInterp=Blue

Din rezultat reise ca avem un obiect GeoTiff, proiectia nu este cunoscutata, dar din coordonatele extend-ului banuim ca este vorba de Stereo70. Acest obiect a fost produs si exportat in suita de softuri ERDAS. Obiectul este un raster color, cu trei benzi RGB, necomprimat, iar sisteml de coordonate din header nu este declarat sau este ne-interpretabil conform topologiei EPSG.

Ultimul exemplu realizat pe un raster de elevatie de la NASA:

$ wget ftp://e0srp01u.ecs.nasa.gov/srtm/version2/SRTM3/Eurasia/N46E038.hgt.zip
$ unzip -e N46E038.hgt.zip
$ gdalinfo N46E038.hgt
Driver: SRTMHGT/SRTMHGT File Format
Files: N46E038.hgt
Size is 1201, 1201
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9108"]],
AXIS["Lat",NORTH],
AXIS["Long",EAST],
AUTHORITY["EPSG","4326"]]
Origin = (71.999583333333334,1.000416666666667)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  (  37.9995833,  47.0004167) ( 37d59'58.50"E, 47d 0'1.50"N)
Lower Left  (  37.9995833,  45.9995833) ( 37d59'58.50"E, 45d59'58.50"N)
Upper Right (  39.0004167,  47.0004167) ( 39d 0'1.50"E, 47d 0'1.50"N)
Lower Right (  39.0004167,  45.9995833) ( 39d 0'1.50"E, 45d59'58.50"N)
Center      (  38.5000000,  46.5000000) ( 38d30'0.00"E, 46d30'0.00"N)
Band 1 Block=1201×1 Type=Int16, ColorInterp=Undefined
NoData Value=-32768
Unit Type: m

Avem un set de date SRTM, format HGT, sistem WGS84, unitate de masura grade, o singura banda pe 16 biti, rezolutia spatiala de 0.00083 grade per pixel.

gdal_translate

Tool pentru conversie inter-formate raster, cu capabilitati de decupaj dupa extend-ul specificat de utilizator, fie in pixeli, fie in unitatile proiectiei obiectului raster dat. Un dezavantaj la acest tool este ca opereaza pe un singur obiect/fisier inputat, lucru logic daca avem in vedere menirea in sine a utilitarului. Sa-l analizam mai de aproape:

$ gdal_translate
Usage: gdal_translate [—help-general]
[-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64 CInt16/CInt32/CFloat32/CFloat64}] [-strict]
[-of format] [-b band] [-outsize xsize] ysize[]
[-scale [src_min src_max [dst_min dst_max]]]
[-srcwin xoff yoff xsize ysize] [-projwin ulx uly lrx lry]
[-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]
[-gcp pixel line easting northing [elevation]]*
[-mo "META-TAG=VALUE"]* [-quiet] [-sds]
[-co "NAME=VALUE"]*
src_dataset dst_dataset

In linii mari, utiliarul face urmatoarele actiuni:

• converteste dintr-un format in altul, acceptand parametri specifici formatului de iesire;
• decupeaza dupa un extend/fereastra specificat de utilizator;
• adauga puncte de control (georeferentiere) pe raster;
• re-asigneaza adancime de biti, chiar prin interpolare daca este necesar;
• asigneaza valori de transparenta;
• manipuleaza benzile ce compun rasterul (le scoate doar daca e necesar);
• manipuleaza interpretarea benzilor tematica/atematica/elevatie/infrarosu ce compun rasterul.

Vom merge mai departe cu exemple concrete. In exemplul al doilea de la gdalinfo am remarcat ca obiectul GeoTiff studiat nu avea informatie spatiala si nu era comprimat. Cu ajutorul unor comenzi simple putem corecta aceste aspecte:

$ gdal_translate -a_srs EPSG:31700 RGB321-L-34-020.tif -co compress=lzw -of GTiff  L-34-020.TIF 
Input file size is 2699, 2659
0…10…20…30…40…50…60…70…80…90…100 %u2013 done.

Prima constatare este legata de viteza si performanta de executie. Din experienta personala am constatat cu stupoare ca obiecte multi Giga-byte pot fi manipulate cu o usurinta si o viteza fantastica, incomparabila cu al unui vendor comercial, ce ofera o interfata mult prea lenesa pentru asa ceva.

Ce s-a obtinut? Sa vedem cu gdalinfo:

$ gdalinfo L-34-020.TIF
Driver: GTiff/GeoTIFF
Files: L-34-020.TIF
Size is 2699, 2659
Coordinate System is:
PROJCS["unnamed",
GEOGCS["Dealul Piscului 1970",
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","31700"]]
Origin = (235438.213641703130000,691205.184220034980000)
Pixel Size = (14.590382532428649,-14.590382532428643)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=IMAGINE TIFF Support
Copyright 1991 – 1999 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: etif.c $ $Revision: 1.10.1.9.1.9.2.11 $ $Date: 2004/09/15 18:42:01EDT $
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
TIFFTAG_RESOLUTIONUNIT=1 (unitless)
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  235438.214,  691205.184) 
Lower Left  (  235438.214,  652409.357) 
Upper Right (  274817.656,  691205.184) 
Lower Right (  274817.656,  652409.357) 
Center      (  255127.935,  671807.271) 
Band 1 Block=2699×1 Type=Byte, ColorInterp=Red
Band 2 Block=2699×1 Type=Byte, ColorInterp=Green
Band 3 Block=2699×1 Type=Byte, ColorInterp=Blue	

Primul avantaj este declararea clara in obiect a proiectiei: EPSG:31700 (Dealul Piscului 1970). Orice copie a fisierului, data oricui, va fi cu siguranta interpretata corect de catre orice aplicatie GIS compatibila EPSG. Pentru mai multe detalii referitoare la proiectiile incluse in baza de date EPSG va invitam sa sa consultati http://www.spatialreference.org, un motor web ce permite cautare prin baza de date, fiind extrem de util pentru clarificarea codului si denumirii unui anumit sistem de coordonate.

Sa vedem putin mai departe diferenta de marime in obiectul rezultat:

$ ls -lh RGB321-L-34-020.tif L-34-020.TIF 
—————+ 1 cbalint        None 17M Jan 30 16:21 L-34-020.TIF
rwx——-+ 1 Administrators None 22M Aug 25 19:59 RGB321-L-34-020.tif

Am comprimat obiectul raster folosind algoritmul LZW (zip) si am castigat 5Mb spatiu de stocare, fara sa alteram vreun pixel, LZW fiind o compresie fara pierdere de calitate (lossless).

Alternativ, am fi putut folosi compresia DEFLATE sau chiar JPEG. Sa vedem ce obtinem cu DEFLATE (o varianta zip mai agresiva):

$ gdal_translate -a_srs EPSG:31700 RGB321-L-34-020.tif -co compress=deflate -co zlevel=9 -of GTiff  L-34-020.TIF
$ ls -lh RGB321-L-34-020.tif L-34-020.TIF 
—————+ 1 cbalint        None 15M Jan 30 16:25 L-34-020.TIF
rwx——-+ 1 Administrators None 22M Aug 25 19:59 RGB321-L-34-020.tif

De data aceasta castigul este de 7Mb, ceea ce este semnificativ, avand in vedere ca am putea avea sute de rasteri de acest gen, sau poate cateva mii. Problema este ca DEFLATE, in comparatie cu LZW, nu este interpretat si acceptat de toate aplicatiile GIS existente. Atentie deci la conversii, GDAL poate face multe minuni dar nu este garantata interoperabilitatea cu alte aplicatii, in speta cu cele proprietare. De aceea este bine de realizat un test prealabil, care sa asigure utilizatorul ca formatul produs este in regula si satisfacator pentru stocare si manipulare. Pentru a vedea mai multi parametri acceptati de driverul GeoTiff, sau alte posibile drivere de iesire, puteti apela manualul inclus sau cel online (ex.: pt geotiff, http://www.gdal.org/frmt_gtiff.html), unde se mentioneza pentru fiecare format compresile acceptate. Pentru GeoTiff mai remarcam CCIT, o compresie folosita in facsimil, potrivita pentru benzile alb/negru, cu adancime 1bit, putandu-se obtine compresii „lossless” la rate extrem de mari.

Un alt exemplu ar fi sa georeferentiem plansa raster din primul exemplu de la sectiunea gdalinfo. Pentru acest lucru vom avea nevoie de o aplicatie de editare a imaginilor, cu ajutorul careia vom afla coordonatele in pixeli ale punctelor de control folosite la georeferentierea hartii. Personal am folosit GIMP, dar orice alt soft de manipulare a imaginilor raster este bun. Probabil, in mod ironic, in Paintbrush sau Irfanview se puteau depista mai rapid coordonatele in pixeli ale rasterului. Pentru acest pas se poate dezvolta o tehnica foarte rapida, care sa ne permita automatizarea proceseului. De exemplu GIMP ar putea fi folosit pentru a deschide un calup de imagini ne-referentiate si adaugarea a cate unui pixel singular, de o anumita culoare, sa zicem in chiar in layerul de transparenta, pixel care cu un tool din command line poate fi reperat si i se poate asocia o coordonata calculata dupa logica de numerotare a planselor. Aici incap multe discutii dar nu vom intram in detalii ci doar vom concluziona ca procesul se poate automatiza.

Revenind la exemplul nostru, am constatat in GIMP ca plansa /images/austrians/p_39_47_100.jpg are urmatoarele puncte de control pixel—>coordonate:

• 338, 85 —> 38.5, 47.5
• 3620, 220 —> 39.5, 47.5
• 3442, 5061 —> ,39.5, 46.5
• 98, 4912 —> 38.5, 46.5

Astfel ca putem trece la georeferentierea pozei:

$ gdal_translate p_39_47_100.jpg -gcp 338 85 38.5 47.5 -gcp 3620 220 39.5 47.5 -gcp 3442 5061 39.5 46.5 -gcp 98 4912 38.5 46.5 -of GTiff -co compress=lzw p_39_47_100.TIF
Input file size is 3729, 5139
0…10…20…30…40…50…60…70…80…90…100 – done.

In felul acesta am obtinut un obiect raster in format GeoTiff cu puncte de control:

$ gdalinfo  p_39_47_100.TIF
Driver: GTiff/GeoTIFF
Files: p_39_47_100.TIF
Size is 3729, 5139
Coordinate System is `’
GCP Projection = 
GCP[  0]: Id=1, Info=(338,85) -> (38.5,47.5,0)
GCP[  1]: Id=2, Info=(3620,220) -> (39.5,47.5,0)
GCP[  2]: Id=3, Info=(3442,5061) -> (39.5,46.5,0)
GCP[  3]: Id=4, Info=(98,4912) -> (38.5,46.5,0)
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 5139.0)
Upper Right ( 3729.0,    0.0)
Lower Right ( 3729.0, 5139.0)
Center      ( 1864.5, 2569.5)
Band 1 Block=3729×1 Type=Byte, ColorInterp=Red
Band 2 Block=3729×1 Type=Byte, ColorInterp=Green
Band 3 Block=3729×1 Type=Byte, ColorInterp=Blue

Remarcam ca acest obiect nu contine un sistem de coordonate sau unitate de masura dar contine punctele de control inserate pentru fiecare punct arondat unei valoari in coordonate geografice. GeoTiff-ul a fost comprimat LZW pentru economie de spatiu.

Vom incerca sa asignam si un sistem de coordonate. Pentru aceasta ne vom uitam la motorul EPSG online. Aici am gasit un echivalent pentru sistemul austro-ungar, denumit Ferro si avand codul EPSG:4805: http://www.spatialreference.org/ref/epsg/4805/. Data fiind lipsa detaliilor asupra plansei este destul de greu sa apreciem daca sistemul identificat este cel mai bun posibil, dar cu siguranta pentru exemplul nostru este suficient.

$ gdal_translate p_39_47_100.jpg -gcp 338 85 38.5 47.5 -gcp 3620 220 39.5 47.5 -gcp 3442 5061 39.5 46.5 -gcp 98 4912 38.5 46.5 -of GTiff -co compress=lzw -a_srs EPSG:4805  p_39_47_100.TIF
Input file size is 3729, 5139
0…10…20…30…40…50…60…70…80…90…100 – done.

$ gdalinfo  p_39_47_100.TIF
Driver: GTiff/GeoTIFF
Files: p_39_47_100.TIF
Size is 3729, 5139
Coordinate System is `’
GCP Projection = GEOGCS[“MGI (Ferro)”,DATUM[“unknown”,SPHEROID[“unretrievable – using WGS84”,6378137,298.257223563]],PRIMEM[“Greenwich”,0],UNIT[“degree”,0.0174532925199433],AUTHORITY[“EPSG”,“4805”]]
GCP[  0]: Id=1, Info=(338,85) -> (38.5,47.5,0)
GCP[  1]: Id=2, Info=(3620,220) -> (39.5,47.5,0)
GCP[  2]: Id=3, Info=(3442,5061) -> (39.5,46.5,0)
GCP[  3]: Id=4, Info=(98,4912) -> (38.5,46.5,0)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 5139.0)
Upper Right ( 3729.0,    0.0)
Lower Right ( 3729.0, 5139.0)
Center      ( 1864.5, 2569.5)
Band 1 Block=3729×1 Type=Byte, ColorInterp=Red
Band 2 Block=3729×1 Type=Byte, ColorInterp=Green
Band 3 Block=3729×1 Type=Byte, ColorInterp=Blue

De aceasta data remarcam ca avem inclusa atat proiectia cat si punctele de control, in concluzie avem un obiect raster integru si transferabil oricui fara prea multe note de explicatii asupra fisierului. Cine este curios poate citi inapoi, utilizand gdalinfo, modalitatea in care el a fost georeferentiat si punctele de control utilizate.

Observatie
Nu toate aplicatiile GIS accepta puncte de control, atentie deci la aceasta modalitate de a stoca obiecte spatiale raster. Partea buna la aceasta metoda de inserare a punctelor de control este ca se obtine un singur fisier, in starea nativa si nealterata a imaginii, fara se deformeaze pixelii, la care se adauga in header informatiile spatiale sub forma unor puncte, lasand mai departe la latitudinea utilizatorilor de a-si alege algoritmul preferat si suita GIS preferata de a aduce/exporta rasterul la forma geo-referentiata finala.

gdalwarp

Utilitar pentru reproiectia datelor spatiale. Sa facem cunostinta mai indeaproape cu capabilitatile acestui utilitar:

$ gdalwarp
Usage: gdalwarp [—help-general] [—formats]
[-s_srs srs_def] [-t_srs srs_def] [-order n] ] [-tps] [-et err_threshold]
[-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]
[-wo “NAME=VALUE”] [-ot Byte/Int16/…] [-wt Byte/Int16]
[-srcnodata “value [value…]”] [-dstnodata “value [value…]”] -dstalpha
[-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
[-of format] [-co “NAME=VALUE”]* srcfile* dstfile
Available resampling methods: near (default), bilinear, cubic, cubicspline, lanczos.

Practic gdalwarp permite:

• reproiectarea dintr-un sistem de proiectie in altul, pentru orice format raster de intrare/iesire;
• decuparea dupa o fereastra predefinita;
• re-esantionarea unui raster la o scara mai redusa;
• manipularea transparentei;
• posibilitatea de selectarea a algoritmului de interpolare si re-esantionare a pixelilor in cadrul oricarui proces de manipulare;

Pentru exemplificare vom incerca sa reproiectam plansa SRTM, utilizata in exemplul anterior, din WGS84 (cod EPSG:4326) in Stereo70 (cod EPSG:31700), iar rezultatul il vom salva in format GeoTiff, comprimat LZW, folosind un algoritm de interpolare a pixelilor cat mai fidel, sa zicem lanczos.

$ gdalwarp -t_srs EPSG:31700 -r lanczos -of GTiff -co compress=LZW N46E038.hgt N46E038-st70.TIF
Creating output file that is 1188P x 1541L.
Processing input file N46E038.hgt.
Using internal nodata values (eg. -32768) for image N46E038.hgt.
0…10…20…30…40…50…60…70…80…90…100 – done.

Mai departe vom reproiecta si plansa din sistemul Ferro, cu puncte de control GCP, generata de noi la sectiunea “gdal_translate”, folosind algoritmul lanczos pentru resampling si proiectia Stereo70. Deoarece avem de aface cu puncte de control in plansa noastra vom lasa la latitudinea gdalwarp sa aleaga gradul trasformarii „affine”, de asociere a pixelilor pe domeniu spatial. Cel mai probabil GDAL va alege gradul 1, pentru transformarea affine cu 4 puncte de control, atatea cate am declarat anterior, pe care o va reproduce impecabil, iar calitatea pixelilor adiacenti si propagarea culorilor spre invecinatati (resampling-ul) va fi una net superioara datorita metodei lanczos alese la re-esantionare pe matricea de pixeli. Din pacate, in acest exemplu, calitatea originala a plansei este relativ slaba, ceea ce micsoreaza sansele de a observa efectelele exacte ale sampler-ului lanczos.

$ gdalwarp -t_srs EPSG:31700 p_39_47_100.TIF -r lanczos -of GTiff -co compress=LZW p_39_47_100-st70.TIF
Creating output file that is 4380P x 5596L.
Processing input file p_39_47_100.TIF.
0…10…20…30…40…50…60…70…80…90…100 – done.

Procesul dureaza ceva mai mult datorita sampler-ului lanczos, de inalta calitate, recomandat in special pentru aerofotograme de mare rezolutie, format RGB pe 24 biti. Plansa rezultata, in proiectie Stereo70, este prezentata in Figura 1.

Inca de la prima vedere se pot ramarca unele inconveniente. Este vorba de cele trei chenare in surplus, unul negru, al doilea alb, plus chenarul inutil al plansei insasi, unde apar gradatiile. Chenarul negru apare in urma trasformarii in Stereo70. Asta deoarece rasterul are forma unei fereastre rectangulare, ne-existand fisiere raster cu forma iregulata. Chenarul alb apare din urma scanarii plansei. Ambele pot constitui un obstacol atunci se ridica problema alipirii cu o plansa adiacenta. Ultimul exemplu isi propune, in patru pasi, sa combata elegant aceste neajunsuri cu ajutorul “gdal_translate” si “gdalwarp”.

• georeferentiere imagine

$ gdal_translate p_39_47_100.jpg -gcp 338 85 38.5 47.5 -gcp 3620 220 39.5 47.5 -gcp 3442 5061 39.5 46.5 -gcp 98 4912 38.5 46.5 -of GTiff -co compress=lzw -a_srs EPSG:4805  p_39_47_100.TIF
Input file size is 3729, 5139
0…10…20…30…40…50…60…70…80…90…100 – done.

• eliminarea chenarului cu gradatii, folosind optiunea -te pentru a specifica fereastra de decupaj in grade

$ gdalwarp -te 38.5 46.5 39.5 47.5  p_39_47_100.TIF  p_39_47_100-noborder.TIF
Creating output file that is 4037P x 4037L.
Processing input file p_39_47_100.TIF.
0…10…20…30…40…50…60…70…80…90…100 – done.

• reproiectare in Stereo70

$ gdalwarp  -t_srs EPSG:31700 p_39_47_100-noborder.TIF p_39_47_100-st70.TIF
Creating output file that is 4006P x 5208L.
Processing input file p_39_47_100-noborder.TIF.
0…10…20…30…40…50…60…70…80…90…100 – done.

• asociere culoare negru (cod culoare = #0) ca fiind transparenta

$ gdal_translate -a_nodata 0  p_39_47_100-st70.TIF -co compress=lzw p_39_47_100-st70-noborder.TIF
Input file size is 4006, 5208
0…10…20…30…40…50…60…70…80…90…100 – done.

Verificam rezultatul:

$ gdalinfo p_39_47_100-st70-noborder.TIF 
Driver: GTiff/GeoTIFF
Files: p_39_47_100-st70-noborder.TIF
 p_39_47_100-st70-noborder.aux
 p_39_47_100-st70-noborder.rrd
Size is 4006, 5208
Coordinate System is:
 PROJCS[“unnamed”,
 GEOGCS[“Dealul Piscului 1970”,
 DATUM[“unknown”,
 SPHEROID[“unretrievable – using WGS84”,6378137,298.257223563]],
 PRIMEM[“Greenwich”,0],
 UNIT[“degree”,0.0174532925199433]],
 UNIT[“metre”,1,
 AUTHORITY[“EPSG”,“9001”]],
 AUTHORITY[“EPSG”,“31700”]]
Origin = (1514274.239675097900000,767631.323753058790000)
Pixel Size = (23.757760795673267,-23.757760795673267)
Metadata:
 AREA_OR_POINT=Area
Image Structure Metadata:
 COMPRESSION=LZW
 INTERLEAVE=PIXEL
Corner Coordinates:
 Upper Left  ( 1514274.240,  767631.324) 
 Lower Left  ( 1514274.240,  643900.906) 
 Upper Right ( 1609447.829,  767631.324) 
 Lower Right ( 1609447.829,  643900.906) 
 Center      ( 1561861.035,  705766.115) 
Band 1 Block=4006×1 Type=Byte, ColorInterp=Red
 NoData Value=0
 Overviews: 1002×1302, 501×651, 251×326, 126×163, 63×82
 Metadata:
 LAYER_TYPE=athematic
Band 2 Block=4006×1 Type=Byte, ColorInterp=Green
 NoData Value=0
 Overviews: 1002×1302, 501×651, 251×326, 126×163, 63×82
 Metadata:
 LAYER_TYPE=athematic
Band 3 Block=4006×1 Type=Byte, ColorInterp=Blue
 NoData Value=0
 Overviews: 1002×1302, 501×651, 251×326, 126×163, 63×82
 Metadata:
 LAYER_TYPE=athematic

Observam sintagma NoData Value=0 pe cele 3 benzi RGB, asta inseamna ca orice soft GIS va interpreta valoarea 0 a culorii ca fiind transparenta, cel putin in aplicatiile ESRI si ERDAS functioneaza, rasterul fiind corect afisat (Figura 2). Este mai comod sa atribuim valoarea NoData in headerul rasterului, decat sa definim o a patra banda pe langa RGB si s-o declaram transparenta, lucru ce se va simti in spatiul de stocare necesar, fisierul crescand in marime datorita benzii alpha suplimentare.

Disclaimer

Acest articol este lipsit de orice drept de proprietate intelectuala restrictiva publicului, este declarat cu drepturi depline de uz public, iar unele materiale prezentate se supun drepturilor de autor specificate pe site-ul de provenienta geo-spatial.org. Suita GDAL este sub licenta Public Domain compatibila cu acest articol, iar marcile inregistrate Leica, ESRI, ERDAS, GlobbalMapper, Microsoft, Google se pot supune unor legi prin simpla folosire a denumirii, in aceasta privinta autorul articolului declara ca nu este afiliat sub nici un fel acestor marci, ele fiind amintite si utilizate ca denumire in scop pur documentativ si exemplificatoriu.