Realizarea unui "super mozaic" raster folosind instrumente open source

1. Introducere

Mozaicarea foilor de hartă sau a scenelor satelitare este o activitate des întîlnită în cadrul proiectelor cu tentă geospațială. Avantajul: un singur fișier este mult mai ușor de manipulat decît o colecție de zeci, sute sau chiar mii de fișiere. Dezavantajul: o dată cu creșterea numărului de scene/foi de hartă dintr-un mozaic crește și dimensiunea fișierului. O soluție de compromis este stocarea mozaicului într-un format de fișier special optimizat pentru o compresie cît mai mare a datelor. Printre aceste formate se numără: MrSID (Multi-resolution Seamless Image Database), ECW (Enhanced Compression Wavelet), JPEG2000. În tutorialul de față ne propunem să obținem un mozaic național din foile de hartă din a treia ridicare austriacă. De-a lungul timpului mai mulți membrii ai comunității geo-spatial.org și-au exprimat interesul pentru un astfel de produs. Pentru distribuție vom stoca fișierul în format ECW. Am ales acest format deoarece știam din experiențele anterioare că oferă un raport foarte bun calitate/compresie și, mai ales, datorită faptului că SDK-ul de compresie este publicat sub licență open source. Pentru manipularea datelor am ales cuplul de librării open source GDAL/OGR. Chiar dacă GDAL/OGR rulează în linie de comandă, utilizatorii mai puțin experimentați nu au motive să fie intimidați. Comenzile în sine sînt destul de simple iar eu îmi iau angajamentul solemn de a documenta foarte bine fiecare pas. Pentru a înțelege mai bine utilitarele GDAL vă recomand tutorialul Manipularea datelor spațiale folosind GDAL scris de Cristian Balint. Platforma pentru care voi ilustra procedura de mozaicare este Windows XP. Utilizatorii de sisteme Linux sau MacOSX pot urma aceeași pași, diferențele apărînd la modul de instalare a SDK-ului de compresie. ERMapper (acum ERDAS) oferă versiuni compilate doar pentru platformele Windows. Utilizatorii ce folosesc alte sisteme de operare vor trebui să-și compileze singuri codul sursă.

Avertisment!
Tutorialul va fi destul de lung și, probabil, plictisitor. Am încercat sa-l scriu cît se poate de detaliat pentru a fi accesibil oricărui tip de utilizator.

2. Date

După cum am scris și mai sus, pentru exemplificare vom folosi drept sursă de date hărțile austriece executate de Imperiul Habsburgic în cadrul celei de a treia ridicări topografice militare. Foile de hartă au fost georeferențiate, în sistem Stereo70, folosind puncte de interes comune, identificate pe hărțile austriece și pe cele topografice recente. Puteți descărca individual aceste foi de hartă de la https://www.geo-spatial.org/vechi/download/harile-austriece-1910-reproiectate-in-stereo70. Hărțile sînt oferite în două formate de fișier: GeoTiff și MrSID. Pentru eficiență am optat pentru datele în format GeoTiff. Asta deoarece GDAL/OGR suportă doar citirea datelor în format MrSID, nu și scrierea, iar unul din pașii intermediari presupune scrierea de informație în fișierul sursă. Automat, ar fi trebuit să convertim fișierele MrSID într-un format de fișier în care GDAL/OGR poate și scrie, iar GeoTiff este un asemenea format. După descărcare și dezarhivare veți obține 40 de fișiere cu extensia .tif și încă 40 cu extensia .tfw, după modelul 38_45.tif & 38_45.tfw. Vă recomand să plasați aceste fișiere într-un folder separat.

Notă
Cei care au descărcat anterior hărțile austriece sînt invitați să o facă din nou. În timp ce scriam acest tutorial am depistat o serie de probleme în fișierele .tif originale și am refăcut complet setul de date.

3. Instrumente

În toate etapele parcurse pentru realizarea mozaicului am folosit numai aplicații libere open source. Lista include:

3.1. GDAL/OGR. L-am instalat folosind suita FWTools. FWTools are avantajul de a include mediul de dezvoltare Python și o serie de pachete (bindings) necesare pentru operația de mozaicare. Cei care optează pentru instalarea separată a librăriei GDAL/OGR vor trebui să instaleze separat Python , GDAL Python, NumPy și vor trebui să configureze manual o serie de variabile de sistem. Pentru a testa instalarea GDAL/OGR vă invit să deschideți FWTools Shell dînd dublu click pe scurtătura de pe desktop sau din Start/Programs/FWTools/FWTools Shell. În fereastra nou deschisă navigați către folderul în care ați salvat imaginile (în cazul meu, comenzile au fost (1) h: pentru a schimba drive-ul și (2) cd H:\Proiecte\geo-spatial.org\Austrian pentru a naviga în folderul cu pricina), apoi scrieți gdalinfo 38_45.tif și dați Enter. Rezultatul ar trebui să fie similar (cu excepția căilor) cu cel prezentat în Figura 1.

Figura 1. Testarea instalării GDAL/OGR folosind comanda gdalinfo.

De aici aflăm o serie de informații despre imaginea în cauză, printre care: dimensiune imagine (3934 x 5282 pixeli), dimensiune celula (23m), coordonate colțuri, număr benzi, sistemul de coordonate etc.

3.2. ECW JPEG2000 Codec SDK. Versiunea compilată poate fi descărcată la adresa http://erdas.com/tabid/84/currentid/1142/default.aspx. În prealabil, este necesară crearea unui cont (gratuit) de utilizator pe site-ul ERDAS. Tot de aici, utilizatorii de sisteme de operare non Windows pot descărca codul sursă al SDK-ului. În continuare instalați aplicația folosind opțiunile implicite. O dată finalizat procesul de instalare copiați fișierele .dll (NCScnet.dll, NCSEcw.dll, NCSEcwC.dll, NCSUtil.dll) din folderul C:\Program Files\Earth Resource Mapping\ECW SDK\redistributable\vc71 în C:\WINDOWS\system32.

3.3. Client GIS desktop. Necesar pentru vizualizarea datelor inițiale și a rezultatului final. Personal am mers pe QGIS. La fel de bine puțeți folosi uDig, OpenJUMP, gvSIG sau o altă aplicație capabilă să citească fișiere raster georeferențiate.

3.4. Componente opționale. Pentru crearea unui fișier de comenzi care să permită procesarea în mod batch a setului am mai avea nevoie de un editor de text (pe platforme Windows vă recomand Notepad++), un mediu de programare (alegerea mea: PHP) sau o aplicație de calcul tabelar (recomand aplicația Calc din suita OpenOffice). Comenzile pot fi scrise si manual, fără a necesita aceste instrumente (poate doar cu excepția editorul de text). Doar că timpul de lucru va crește sensibil.

4. Preprocesare

Pentru a fi gata de mozaicare fișierele de intrare au nevoie de o procesare intermediară. Scopul acestui pas este eliminarea informației marginale de pe hartă și păstrarea informației utile. Datorita convențiilor de împărțiere a foilor de hartă, în marea majoritate a cazurilor, în sistemele de coordonate carteziene, suprafață utilă a hărților topografice are forma unui trapez. De aici și denumirea de "trapeze", folosită adesea pentru a desemna foile de hartă topografică. Atunci cînd harta este scanată, apoi georeferențiată în sistemul de coordonate cartezian, iar informația marginală este eliminată vom obține un astfel de trapez. Însă, la stocarea informației (trapezului) într-un format de fișier digital (.bmp, .jpg, .tif, etc.), trapezul va fi încadrat într-un dreptunghi și salvat ca atare pe disc. Aceasta deoarece toate aceste formate de fișiere stochează informația sub formă matricială. În funcție de formatul de fișier și de aplicația folosită, pixelilor rezultați din diferența dintre dreptunghi și trapez li se va aloca o anumită valoare. Pentru rezultate optime la mozaicare, acești pixeli trebuiesc ignorați (asociați cu nodata). Figura 2 prezintă alăturarea hărților cu tot cu informația marginală.

Figura 2. Exemplu de alăturare a două foi de hartă fără eliminarea informație marginale.

Este evident că pentru a obține un produs coerent este necesar să eliminăm chenarul alb și elementele grafice din afara chenarului interior al hărții (Figura 3).

Figura 3. Exemplu de alăturare a două foi de hartă după eliminarea informație marginale.

Deoarece GDAL, librăria cu care ne propunem să preprocesăm datele, nu permite extragerea unui subset dintr-un raster cu o formă geometrică, altă decît cea de dreptunghi cu laturile paralele cu axele sistemului de coordonate, va trebui să apelăm la unele trucuri. Problema dată poate fi rezolvată în două moduri:

4.1. Reproiectînd hărțile în coordonate geografice

Formă de trapez, amintită anterior, este valabilă atîta vreme cît harta este "ținută" în sistemul de coordonate cartezian (Figura 4). Dacă harta este reproiectată într-un sistem de coordonate geografic, forma suprafeței utile din hartă se transformă într-un dreptunghi cu laturile paralele cu rețeaua de meridiane și paralele geografice. De aici înainte putem folosi utilitarul gdalwarp, cu parametrul -te pentru a stabili extinderea maximă a fișierului rezultat.

Figura 4. Exemplu de foaie de hartă austriacă georeferențiată în Stereo70.

Înainte de a rula comanda gdalwarp pentru reproiectare și eliminarea informației marginale avem nevoie să identificăm coordonatele colțurilor suprafeței utile din hartă. Fiecare planșă din seria austriacă are dimensiunea de 1 grad longitudine x 1 grad latitudine. Inspectînd planșa intitulată 38_45, într-un client desktop, vom observa că aceasta se întinde între 37º30'00'' - 38º30'00'' pe longitudine și 44º30'00'' - 35º30'00'' pe latitudine. Valorile longitudinale sînt mai mari decît cele întîlnite în mod obișnuit pe teritoriul României deoarece aceste hărți folosesc drept referință meridianul Ferro și nu Greenwich. Diferența dintre Ferro și Greenwich, calculată pentru hărțile din a treia ridicare topografică austriacă, este de 17º39'46.02'' (Gábor TIMÁR). Pentru a le putea introduce în GDAL, coordonatele au fost convertite din formatul GMS (Grad:Minut:Secundă) în format GZ (Grade Zecimale), obținînd 37.50 - 38.50 cu 44.50 - 45.50.

La reproiectare vom folosi sistemul de coordonate cu codul EPSG 4805. Același sistem a fost folosit și de Cristian Balint în tutorialul său, ilustrat tot cu foile de hartă austriacă.

Pentru realizarea operațiunii vă invit să deschideți o instanță FWTools Shell, să navigați către folderul unde se găsesc hărțile și să scrieți următoarea linie de cod:

	gdalwarp -of GTiff -s_srs "EPSG:31700" -t_srs "EPSG:4805" -te 37.50 44.50 38.50 45.50 38_45.tif 38_46_aus.tif	

În fereastra FWTools Shell, comanda ar trebui să producă următorul rezultat:

	H:\Proiecte\geo-spatial.org\Austrian>gdalwarp -of GTiff -s_srs "EPSG:31700" -t_srs "EPSG:4805" -te 37.50 44.50 38.50 45.50 38_45.tif 38_46_aus.tif
	Creating output file that is 4075P x 4075L.
	Processing input file 38_45.tif.
	0...10...20...30...40...50...60...70...80...90...100 - done.	

Fișierul rezultat ar trebui să fi similar cu cel prezentat în Figura 5.

Figura 5. Rezultatul reproiectării unei foi de hartă în sistem de coordonate geografic și eliminarea informației marginale.

Comanda gdalwarp din GDAL permite utilizatorilor conversia imaginilor dintr-un sistem de coordonate în altul, precum și aplicarea mai multor tipuri de transformări. În cazul nostru, parametrii utilizați împreună cu această comandă au următorul rol:

• -of GTiff specifică formatul de fișier în care va fi stocat rezultatul prelucrării, în cazul nostru GeoTiff. Am fi putut ignora acest parametru deoarece, în mod implicit, GDAL folosește formatul GeoTiff dacă parametrul _-of_ nu a fost declarat. O listă completă a formatelor suportate de GDAL se găsește aici.
• -s_srs 'EPSG:31700' specifică sistemul de coordonate al datelor de intrare (source spatial reference set), în cazul nostru Stereo70, reprezentat prin codul EPSG 31700. Și acest parametru putea fi ignorat deoarece GDAL ar fi găsit această informație în headerul fișierului de intrare. Precizarea este însă foarte importantă atunci cînd avem de-a face cu fișiere în format .tif însoțite doar de fișiere .tfw, fără header de tip GeoTiff sau fișier .prj, sau cu alte tipuri de fișiere ce nu conțin în mod explicit definiția sistemului de coordonate.
• -t_srs 'EPSG:4805' specifică sistemul de coordonate pentru fișierul rezultat din prelucrare (target spatial reference set). După cum am precizat anterior, vom folosi sistemul EPSG 4805, numit și MGI (Ferro).
• -te 37.50 44.50 38.50 45.50 precizează extinderea spațială a fișierului rezultat (format minX minY maxX maxY), folosind valori în sistemul de coordonate țintă.
• 38_45.tif numele fișierului de intrare.
• 38_45_aus.tif numele fișierului de ieșire (fișierul în care va fi stocat rezultatul transformării).

Pentru lista completă, însoțită de descrieri și exemple, a parametrilor ce pot însoți comanda gdalwarp vă recomand pagina http://gdal.org/gdalwarp.html.

În continuare putem aplica același procedeu pentru toate foile de hartă. La sfîrșit, fișierele rezultate vor fi mozaicate iar mozaicul va fi reproiectat din nou în Stereo70 și convertit în format ECW.

Setul de date folosit pentru ilustrare vine cu unele probleme. Fiind vorba de hărți istorice, pentru care nu s-au cunoscut precis parametrii proiecței, o serie de deformări au fost induse în procesul de georeferențire. Procedura de georeferențiere s-a bazat pe identificarea de puncte comune (de obicei biserici și mănăstiri) pe hărțile austriece și pe hărțile topografice militare recente. Precizia și densitatea acestor puncte comune diferă de la o foaie de hartă la alta și de la o regiune istorică la alta. Din această cauză, la unele foi de hartă, precizia georeferențierii este mai slabă, iar anumite puncte, cum ar fi colțurile hărții, nu se mai regăsesc la coordonatele indicate de valorile înscrise pe marginea hărții. Astfel de planșe se găsesc cu predilecție în sudul Munteniei și în Dobrogea. În aceste cazuri, la operațiunea de "tăiere" pe coordonate riscăm să includem mici suprafețe din cadrul exterior al hărții, sau să omitem mici suprafețe din hartă. O ilustrare a acestui gen de problemă este prezentată în Figura 6. Astfel de probleme nu ar fi apărut dacă am fi prelucrat o hartă topografică modernă.

Figura 6. Rezultatul obținut cu gdalwarp pentru planșa 47_45. Liniile roșii unesc colțurile reale (cruce albastră) ale foii de hartă.

4.2. Folosind o mască vectorială

Afirmam mai sus că GDAL "nu știe" să decupeze un raster folosind o mască cu geometrie complexă. Există însă un utilitar, gîndit pentru un cu totul alt scop, care poate fi "convins" să ofere o astfel de funcționalitate. Utilitarul cu pricina se numește gdal_rasterize. Scopul acestuia este de a aplica o mască vectorială peste un raster și de a suprascrie (burn) valorile din celulele "acoperite" de vector cu o valoare nouă. Un exemplu de caz în care o astfel de comandă ar fi utilă este clasificarea unei imagini satelitare acoperită parțial de nori. Într-un astfel de caz putem realiza o mască vectorială pentru zonele acoperite cu nori și pe baza acesteia, folosind comanda gdal_rasterize, să marcăm acele zone din imagine ca nodata. Comanda poate fi rulată folosind mai mulți parametri. Lista completă se găsește aici. Unul din parametrii care face comanda interesantă pentru cazul nostru este -i (invert). Acesta inversează efectul funcției și modifică toate valorile celulelor care nu sînt acoperite de masca vectorială. Scenariul pentru cazul nostru are următorii pași:

• Crearea unei măști vectoriale, de tip poligon, a suprafeței utile din hartă;
• Aplicarea comenzii gdal_rasterize pentru celulele din afara măștii și înlocuirea valorilor existente cu o singură valoare;
• declararea acelei valori ca fiind nodata.

Această abordare ne permite să lucrăm direct în Stereo70.

4.2.1. Crearea măștii se poate face prin vectorizarea manuală a suprafeței utile din hartă (Figurile 7-8) sau prin generare automată, dacă se cunosc principiile modului de împărțire a foilor de hartă.

Figura 7. Vectorizarea suprafeței utile dintr-o planșă.

Figura 8. Vectorizarea suprafeței utile dintr-o planșă (detaliu).

În cazul unor hărți pentru care se pot obține automat coordonatele colțurilor vă prezentăm o metodă rapidă de a le transforma în format ESRI Shapefile pentru a putea fi folosite împreună cu comanda gdal_rasterize. Voi face exemplificarea plecînd de la coordonatele colțurilor (stînga-sus, dreapta-sus, dreapta-jos, stînga-jos) foii de hartă 38_45:

	96093.312685, 457378.714379
	174575.972863, 452947.029426
	169487.214594, 339986.331632
	89574.9434037, 344549.039998	

Una din cele mai simple metode de a transforma aceste coordonate într-un poligon este crearea unui fișier de tip Atlas BNA (.bna). Acesta este un format ASCII cu o structură extrem de simplă, similară cu cea prezentată mai jos:

	"atribut 1","atribut 2",număr vertecși
	x1,y1
	x2,y2
	x3,y3
	...
	
	"atribut 1","atribut 2",număr vertecși
	x1,y1
	x2,y2
	x3,y3
	...

Pentru fiecare element vom avea un header de tipul "atribut 1","atribut 2",număr vertecși, urmat de coordonatele vertecșilor. Respectînd cele scrise mai sus, fișierul .bna pentru planșa 38_45 va arăta așa:

	"38_45", "38_45", 5
	96093.312685, 457378.714379
	174575.972863, 452947.029426
	169487.214594, 339986.331632
	89574.9434037, 344549.039998
	96093.312685, 457378.714379

Observați că prima pereche de coordonate se repetă la finalul fișierului. Acest lucru permite interpretarea elementului ca fiind de tip poligonal. Valorile celor două cîmpuri de atribute nu sînt importante. Am scris ceva acolo pentru a respecta formatul fișierului. Deoarece GDAL/OGR nu permite utilizarea directă măștilor în format .bna în comanda gdal_rasterize, am fost nevoit să convertesc fișierul .bna în formatul de fișier ESRI Shapefile. Acest lucru se poate face foarte simplu cu ajutorul utilitarului ogr2ogr din suita GDAL/OGR:

ogr2ogr -f "ESRI Shapefile" -s_srs "EPSG:31700" -t_srs "EPSG:31700" 38_45.shp 38_45.bna

4.2.2. Avînd la dispoziție masca vectorială, putem trece mai departe la aplicarea efectivă a comenzii gdal_rasterize:

gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 38_45 38_45.shp 38_45.tif

În fereastra FWTools Shell, comanda ar trebui să producă următorul rezultat:

H:\Proiecte\geo-spatial.org\Austrian\de_uploadat>gdal_rasterize -b 1 -b 2 -b 3 -
burn 111 -burn 111 -burn 111 -i -l 38_45 38_45.shp 38_45.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

Atenție!
gdal_rasterize va aplica modificările direct în fișierul de intrare. Dacă doriți să-l păstrați și în formatul original va trebui să faceți o copie de siguranță.

Rezultatul asupra fișierului de intrare este prezentat în Figura 9. După cum se poate observa, toți pixelii din zona exterioară a hărții au acum culoarea gri.

Figura 9. Rezultatul aplicării gdal_rasterize.

În continuare vom lua pe rînd parametrii cu care am rulat comanda:

• -b 1 -b 2 -b 3 Parametrul -b permite selectarea benzii pe care vrem să o modificăm. Prin introducerea repetată a parametrului putem selecta un număr mai mare de benzi. În cazul nostru, fișierul de intrare este compus din trei benzi cu valori de tip byte (între 0 și 255), fiecare bandă corespunzînd uneia din componentele paletei RGB (Red = roșu, Green = verde, Blue = albastru).
• burn 111 -burn 111 -burn 111 Permite introducerea valorii care se va scrie, pentru fiecare pixel, în fiecare bandă în parte. Am ales să scriem valoarea 111 pe fiecare bandă a imaginii. Compuse în format RGB acestea vor da acea nuanță de gri din fundal. Întrebarea care se pune în continuare este: de ce am ales să scriem tocmai valoarea 111? Răspunsul scurt ar fi: pentru că este o valoare care se răgăsește rar în interiorul hărții. La pasul următor, unde vă voi arata cum vom proceda pentru a asocia pixelii din exteriorul hărții cu valoarea nodata, vom vedea de ce 111 este o valoare potrivită pentru acest set de date.
• -i Inversează comportamentul comenzii și aplică modificările în afara măștii vectoriale.
• -l 38_45 Indică stratul care se va folosi ca mască. În cazul formatului ESRI Shapefile, fiecare fișier conține un singur strat. Există însă și formate de tip colecție, care conțin mai multe straturi de informație, iar parametrul -l ne permite sa-l selectăm pe cel care conține masca.
• 38_45.shp Fișierul ce conține stratul mască.
• 38_45.tif Fișierul raster în care se vor aplica modificările.

4.2.2. Urmează să definim prixelii care vor fi ignorați la mozaicare (cei din exteriorul suprafeței utile a hărții). Și aici avem două posibilități:

gdal_translate -a_nodata 111 -co tfw=yes 38_45.tif 38_45_noborder.tif

Parametrii comenzii sînt:

• -a_nodata 111 specifică valoarea care va fi tratată drept nodata.
• -co tfw=yes un nou parametru opțional ce permite crearea unui fișier de tip World File (.tfw). Este util pentru aplicațiile care nu știu să citească informațiile geospațiale din headerul GeoTiff.
• 38_45.tif numele fișierului de intrare.
• 38_45_noborder.tif numele fișierului de ieșire.

Aceasta metodă are un și un dezavantaj. După cum am mai spus, fiecare bandă conține valori cuprinse între 0 și 255. Automat, valorile pe care noi le putem scrie în raster, cu ajutorul gdal_translate, trebuie sa fie și ele cuprinse tot între 0 și 255. Pînă acum nimic grav. Problemele apar la asocierea valorilor 111 cu nodata. Aceasta asociere nu va ține cont de faptul ca pixelii cu valoarea 111 se găsesc în afara hărții sau în interior. Din această cauză, o parte din pixelii valizi din interiorul hărții, care au pe una din cele trei benzi asociată valoarea 111, vor deveni și ei nodata. Comanda standard de definire a unei valori ca fiind nodata face acest lucru pe fiecare banda in parte și nu permite asocieri mai complexe, care să țină cont de valorile unui pixel în toate cele trei benzi. De aceea am căutat o valoare care se regăsește cît mai rar în pixelii din interiorul hărții.

A doua metodă presupune crearea unei noi benzi, de tip alpha. Metoda are avantajul că permite selectarea mai precisă a pixelilor nodata, țînînd cont de valorile din toate cele trei benzi. Astfel, utilizatorul poate introduce o combinație de valori RGB și numai pixelii care respectă strict acea combinație vor fi definiți ca nodata. După cum am scris mai sus, metoda aceasta va crea o nouă bandă. În această bandă vom avea două valori: 0 (pentru pixelii transparenți) și 255 (pentru pixelii opaci). Dezavantajul acestei metode este că mărește dimensiunea fișierului. De asemnea, nu toate aplicațiile GIS de pe piață știu să interpreteze banda de transparență. Mai lentă și mare consumatoare de spațiu, metoda se adresează celor "chițibușari", interesați de o calitate cît mai bună a produsului final. Comanda care produce acest rezultat are forma de mai jos:

gdalwarp -dstalpha -srcnodata "111 111 111" -dstnodata 0 -co tfw=yes 38_45.tif 38_45_noborder.tif

Efectul: doar pixelii ce respectă combinația R=111; G=111; B=111 vor fi asociați cu nodata. Explicația parametrilor:

• -dstalpha creează o nouă bandă de tip alpha.
• -srcnodata "111 111 111" specifică combinația RGB pentru pixelii nodata.
• -dstnodata 0 scrie valoarea 0 în banda alpha pentru toți pixelii identificați ca R=111;B=111;G=111. După cum am precizat anterior, valoarea 0 în banda alpha identifică pixelii transparenți.
• -co tfw=yes aceeași explicație ca la comanda anterioară.
• 38_45.tif numele fișierului de intrare.
• 38_45_noborder.tif numele fișierului de ieșire.

Rezultatul definirii pixelilor nodata, indiferent de metoda aleasă ar trebui să producă, la vizualizarea într-un client GIS, efectul ilustrat în Figura 10.

Figura 10. Rezultatul aplicării transparenței.

Procedura se repetă pentru restul foilor de hartă.

La fel ca și la metoda bazată pe reproiectarea hărților într-un sistem de coordonate geografic, particularitățile amintite ale setului de date utilizat pentru ilustrare, determină unele anomalii pe foile de hartă din sudul țării. Figura 11 ilustrează acest lucru tot pentru planșa 47_45, din Delta Dunării.

Figura 11. Aplicarea măștii vectoriale (ilustrată cu roșu) peste planșa 47_45.

5. Mozaicare

O dată preprocesate, hărțile noastre sînt gata de mozaicare. Nu cred că veți fi mirați dacă vă spun că și pentru acest proces GDAL oferă două metode de lucru.

5.1. gdal_merge.py

gdal_merge.py este un script, scris în Python, care permite mozaicarea rapidă a fișierelor raster. Dezavantajul major al comenzii este că "nu știe" să interpreteze banda alpha, de transparență. Cu ajutorul parametrului -n se poate defini o valoare nodata la nivel de bandă. O selecție complexă, folosind valorile unui pixel în toate cele trei benzi, se poate realiza doar prin programare. gdal_merge.py este recomandat pentru cei care au ales să-și preproceseze datele folosind gdal_translate. Pentru setul nostru de date, comanda va avea următoarea formă:

gdal_merge -of Gtiff -o merge.tif -co tfw=yes -co bigtiff=yes -n 111 38_45_noborder.tif 38_46_noborder.tif 39_45_noborder.tif 39_46_noborder.tif 39_47_noborder.tif 40_44_noborder.tif 40_45_noborder.tif 40_46_noborder.tif 40_47_noborder.tif 40_48_noborder.tif 41_44_noborder.tif 41_45_noborder.tif 41_46_noborder.tif 41_47_noborder.tif 41_48_noborder.tif 42_44_noborder.tif 42_45_noborder.tif 42_46_noborder.tif 42_47_noborder.tif 42_48_noborder.tif 43_44_noborder.tif 43_45_noborder.tif 43_46_noborder.tif 43_47_noborder.tif 43_48_noborder.tif 44_44_noborder.tif 44_45_noborder.tif 44_46_noborder.tif 44_47_noborder.tif 44_48_noborder.tif 45_44_noborder.tif 45_45_noborder.tif 45_46_noborder.tif 45_47_noborder.tif 45_48_noborder.tif 46_44_noborder.tif 46_45_noborder.tif 46_46_noborder.tif 46_47_noborder.tif 47_45_noborder.tif

Explicația parametrilor:

• -of Gtiff specifică formatul pentru fișierul nou rezultat. Dacă parametrul nu este specificat, gdal_merge.py va folosi în mod implicit formatul GeoTiff (am ales sa-l includ totuși în comandă pentru o mai bună înțelegere).
• -o merge.tif numele fișierului de ieșire.
• -co tfw=yes creează fișier .tfw.
• -n 111 valoarea pixelilor nodata, ignorați la mozaicare.
• 38_45_noborder.tif 38_46_noborder.tif ... lista fișierelor de intrare.

În funcție de puterea mașinii (în special procesorul), procesul poate dura între 15 minute și cîteva ore. Rezultatul: un fișier mozaic, de aproximativ 2.5GB (34778 x 24526 px), ce conține toate foile de hartă topografice austriece. Vizualizat în QGIS, mozaicul arată ca în Figura 12.

Figura 12. Rezultatul mozaicării fișierelor.

5.2. gdalwarp

gdalwarp poate fi folosit și el pentru mozaicare. Comanda arată cam așa:

gdalwarp -of Gtiff -co tfw=yes -co 38_45_noborder.tif 38_46_noborder.tif 39_45_noborder.tif 39_46_noborder.tif 39_47_noborder.tif 40_44_noborder.tif 40_45_noborder.tif 40_46_noborder.tif 40_47_noborder.tif 40_48_noborder.tif 41_44_noborder.tif 41_45_noborder.tif 41_46_noborder.tif 41_47_noborder.tif 41_48_noborder.tif 42_44_noborder.tif 42_45_noborder.tif 42_46_noborder.tif 42_47_noborder.tif 42_48_noborder.tif 43_44_noborder.tif 43_45_noborder.tif 43_46_noborder.tif 43_47_noborder.tif 43_48_noborder.tif 44_44_noborder.tif 44_45_noborder.tif 44_46_noborder.tif 44_47_noborder.tif 44_48_noborder.tif 45_44_noborder.tif 45_45_noborder.tif 45_46_noborder.tif 45_47_noborder.tif 45_48_noborder.tif 46_44_noborder.tif 46_45_noborder.tif 46_46_noborder.tif 46_47_noborder.tif 47_45_noborder.tif merge.tif

Observație
Numele fișierului de ieșire se va indica ultimul.

Față de metoda prezentată anterior are avantajul de a interpreta corect banda de transparență. Din păcate este extrem de lent atunci cînd avem de mozaicat un număr mare de fișiere. În cazul datelor noastre am estimat că ar fi nevoie de aproximaiv 50 de ore pentru a realiza mozaicul cu această comandă (n-am avut rabdare/timp sa-l las pînă la capăt). Recomandarea mea este să folosiți gdalwarp doar atunci cînd aveți de mozaicat un număr mic de fișiere (indiferent de dimensiunea acestora). Altfel, nu prea merită.

Atenție!
Pentru mozaicurile foarte mari (peste 4GB), indiferent de metoda de mozaicare aleasă, este necesară alegerea unui format de fișier fără limită de de dimensiune (cum sînt GeoTiff sau ERDAS Img). Lista completă a formatelor suportate de GDAL, însoțite de restricțiile de dimensiune, poate fi consultată aici. Dacă optați pentru GeoTiff va trebui să folosiți parametrul -co bigtiff=yes. BigTiff este o variantă a formatului Tiff care permite stocarea a mai mult de 4GB de informație. În mod normal, GDAL estimează dimensiunea finală a mozaicului și activează singur aceast parametru. Atunci cînd se folosesc algoritmi de compresie, estimarea GDAL poate fi greșită și este mult mai "sănătos" să specificați manual faptul că fișierul de ieșire este de tip BigTiff.

6. Optimizare

6.1. Pre-generarea de versiuni reeșantionate

Dacă ați replicat cu succes instrucțiunile prezentate pînă în acest punct, ați putut observa că simpla încărcare a mozaicului în QGIS se face în circa 10-20 de secunde. Acest lucru se datoreaza faptului că QGIS este nevoit să creeze o imagine micșorată a fișierului de 2.5GB, pentru a o afișa pe ecran. În acest caz, cel mai întîlnit procedeu de optimizare este pre-generarea de versiuni redimensionate (pyramids) ale fișierului. Numărul de nivele pentru care se creează astfel de piramide depinde de dimensiunile fișierului de intrare. Clientul GIS va ști să foloseească aceste versiuni redimensionate atunci cînd fișierul este vizualizat la nivele de zoom mici, crescînd exponențial viteza de randare. Comanda GDAL cu care putem crea piramide este gdaladdo:

gdaladdo -r gauss -ro merge.tif 2 4 8 16 32 64 128

Parametrii folosiți indică:

• -r gauss metoda de reeșantionare folosită. Alte metode disponibile sînt: nearest, average, cubic, average_mp, average_magphase, mode.
• -ro specifică faptul că piramidele vor fi create într-un fișier extern. Formatul GeoTiff permite stocarea piramidelor și în fișierul principal (.tif).
• merge.tif fișierul de intrare.
• 2 4 8 16 32 64 128 nivelele pentru care se vor crea versiuni reeșantionate. Valorile semnifică factorul cu care se va reduce dimensiunea imaginii. Astfel, prima piramidă este de două ori mai mică decît originalul, următoarea de 4 ori și tot așa pînă la de 128 de ori mai mică.

În urma rulării acestei comenzi s-a obținut un fișier nou, numit merge.tif.ovr, care conține versiunile reeșantionate ale mozaicului. Reîncărcați mozaicul în QGIS și veți observa viteza net superioară cu care fișierul este afișat la nivele mici de zoom. De asemenea, se poate observa, că la aceleași nivele mici de zoom, calitatea imaginii este mult mai bună (Figura 13). Acest lucru se datorează folosirii metodei gauss la reeșantionare. În lipsa piramidelor, QGIS ar fi folosit un algoritm gen nearest neighbor, cu rezultate vizuale mult mai modeste.

Notă
Acest tip de optimizare este util pentru cei care vor prefera să utilizeze mozaicul în format GeoTiff. Cei care vor merge mai departe și vor face conversia în format ECW nu au nevoie de acest lucru. La generarea fișierului .ecw, intern, vor fi create și fișierele piramidă.

Figura 13. Vizualizarea mozaicului în QGIS după generarea piramidelor.

6.2. Compresia datelor în format GeoTiff

Dimensiunea fișierului GeoTiff poate fi redusă și fără a apela la alte formate de fișier. Astfel, intern, informațiile stocate în fișierul GeoTiff pot fi compresate folosind diverși algoritmi. De exemplu, comanda:

gdal_translate -co compress=lzw merge.tif merge_lzw.tif

va micșora dimensiunea fișierului de la aprox. 2.5GB la aprox. 1.7GB, folosind algoritmul de compresie LZW (Lempel–Ziv–Welch).

Atenție!
Nu aplicați o astfel de compresie cînd efectuați operațiunea de mozaicare. Dintr-un motiv inexplicabil, fișierul generat de comanda gdal_merge.py + parametrul -co compress=lzw, a avut dimensiunea de 24GB. Același lucru s-a întîmplat cînd am aplicat și un alt algoritm de compresie (deflate). Voi sesiza această problemă pe lista de discuții GDAL.

6.3. Eliminarea fundalului negru

O altă acțiune de îmbunătățire a mozaicului, cu rol mai mult estetic, este eliminarea fundalului negru rezultat după operațiunea de mozaicare. În mod implicit, GDAL va "umple" acești pixeli marginali cu valoarea 0, pe fiecare bandă RGB. Este suficient să declarăm valorea zero ca fiind nodata pentru a elimina negrul acesta deranjant la ochi:

gdal_translate -a_nodata 0 merge.tif merge_noborder.tif

Rezultatul ar trebui să fie similar cu cel prezentat în Figura 14

Figura 14. Eliminarea fundalului negru din mozaic.

7. Compresie ECW

Conversia în format ECW se face cu ajutorul gdal_translate:

gdal_translate -of ECW -co "LARGE_OK=YES" merge.tif merge.ecw

Datorită unor restricții legate de modul de licențiere a SDK-ului ECW și a librăriei GDAL, în mod implicit, GDAL nu va compresa fișierele mai mari de 500MB. Pentru a înlătura aceasta restricție am folosit parametrul -co "LARGE_OK=YES".

Fișierul astfel obținut este de mai bine de 5 ori mai "condensat" față de varianta în format GeoTiff. În mod implicit, nivelul de compresie este de 75%. Valoarea poate fi crescută sau scăzută folosind parametrul -co target=x, unde x este noua reprezintă noul procent de compresie.

8. Automatizare

Aplicarea manuală a acestor comenzi pentru toate cele 40 de fișiere de intrare vă va "mînca" ceva timp. În continuare vă voi prezenta cîteva metode cu ajutorul cărora, utilizatorii care nu știu cum să folosească clasele și metodele GDAL în medii de dezvoltare ca C++ sau Python, vor putea să genereze rapid un fișier de tip Windows Batch.

Pentru început vom ilustra seria de operațiuni ce trebuiesc făcute pentru preprocesarea fiecărei foi de hartă în parte:

• aplicarea gdal_rasterize pentru asocierea pixelilor ce cad in afara măștii vectoriale cu valoarea 111;
• aplicarea gdal_translate pentru declararea valorii 111 ca nodata;

Pentru prima foaie de hartă comenzile arată așa:

gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 38_45 38_45.shp 38_45.tif
gdal_translate -a_nodata 111 -co tfw=yes 38_45.tif 38_45_noborder.tif

8.1. Metoda 1 (nimic automat)

Deschideți un editor de text, copiați aceste două linii de 40 de ori. Folosiți copy&paste pentru a înlocui "38_45" cu numele celorlalte fișiere. Salvați fișierul ca "procesare_harti.bat". Procedura poate fi aplicată cu succes în cazul datelor noastre. Mai greu ar fi pentru un set de date cu sute sau mii de fișiere de intrare.

8.2. Metoda 2 (semi-automată)

Deschideți o aplicație de calcul tabelar (Ex: OpenOffice Calc, Microsoft Office Excel). Pe coloana A scrieți, pe o linie separată, numele fiecărui fișier (fără extensie). Numele le puteți obține automat cu ajutorul Total Commander (selectați fișierele/Mark/Copy Selected Names to Clipboard) sau din lista de mai jos:

	38_45
	38_46
	39_45
	39_46
	39_47
	40_44
	40_45
	40_46
	40_47
	40_48
	41_44
	41_45
	41_46
	41_47
	41_48
	42_44
	42_45
	42_46
	42_47
	42_48
	43_44
	43_45
	43_46
	43_47
	43_48
	44_44
	44_45
	44_46
	44_47
	44_48
	45_44
	45_45
	45_46
	45_47
	45_48
	46_44
	46_45
	46_46
	46_47
	47_45

Selectați celula B1, plasați cursorul în bara Input line și scrieți comanda:

=CONCATENATE("gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l ";A1;" ";A1;".shp";" ";A1;".tif")

Dați Enter. În celula B1 ar trebui să aveți valoarea:

gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 38_45 38_45.shp 38_45.tif

Selectați din nou celula B1, plasați cursorul în dreptul punctului negru din partea dreaptă-jos a celulei și trageți de acesta pînă ajungeți la linia 40. Daca totul a funcționat corect ar trebui să obțineți valori similare cu cele prezentate în Figura 15.

Figura 15. Crearea unui fișier Windows Batch folosind OpenOffice Calc.

Copiați valorile din primele 40 de linii ale coloanei B în clipboard și transferați-le într-un editor de text. Salvați fișierul ca "procesare_harti.bat".

Repetați operațiunea, de data asta calculînd valorile din coloana B folosind relația:

=CONCATENATE("gdal_translate -a_nodata 111 -co tfw=yes ";A1;".tif";" ";A1;"_noborder.tif")

Copiați noile valori din coloana B și adăugațile la sfîrșitul fișierului "procesare_harti.bat". Dacă totul a mers bine, fișierul trebuie să aibă următorul conținut:

gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 38_45 38_45.shp 38_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 38_46 38_46.shp 38_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 39_45 39_45.shp 39_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 39_46 39_46.shp 39_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 39_47 39_47.shp 39_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 40_44 40_44.shp 40_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 40_45 40_45.shp 40_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 40_46 40_46.shp 40_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 40_47 40_47.shp 40_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 40_48 40_48.shp 40_48.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 41_44 41_44.shp 41_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 41_45 41_45.shp 41_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 41_46 41_46.shp 41_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 41_47 41_47.shp 41_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 41_48 41_48.shp 41_48.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 42_44 42_44.shp 42_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 42_45 42_45.shp 42_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 42_46 42_46.shp 42_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 42_47 42_47.shp 42_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 42_48 42_48.shp 42_48.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 43_44 43_44.shp 43_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 43_45 43_45.shp 43_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 43_46 43_46.shp 43_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 43_47 43_47.shp 43_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 43_48 43_48.shp 43_48.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 44_44 44_44.shp 44_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 44_45 44_45.shp 44_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 44_46 44_46.shp 44_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 44_47 44_47.shp 44_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 44_48 44_48.shp 44_48.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 45_44 45_44.shp 45_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 45_45 45_45.shp 45_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 45_46 45_46.shp 45_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 45_47 45_47.shp 45_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 45_48 45_48.shp 45_48.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 46_44 46_44.shp 46_44.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 46_45 46_45.shp 46_45.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 46_46 46_46.shp 46_46.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 46_47 46_47.shp 46_47.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l 47_45 47_45.shp 47_45.tif
gdal_translate -a_nodata 111 -co tfw=yes 38_45.tif 38_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 38_46.tif 38_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 39_45.tif 39_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 39_46.tif 39_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 39_47.tif 39_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 40_44.tif 40_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 40_45.tif 40_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 40_46.tif 40_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 40_47.tif 40_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 40_48.tif 40_48_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 41_44.tif 41_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 41_45.tif 41_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 41_46.tif 41_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 41_47.tif 41_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 41_48.tif 41_48_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 42_44.tif 42_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 42_45.tif 42_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 42_46.tif 42_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 42_47.tif 42_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 42_48.tif 42_48_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 43_44.tif 43_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 43_45.tif 43_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 43_46.tif 43_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 43_47.tif 43_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 43_48.tif 43_48_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 44_44.tif 44_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 44_45.tif 44_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 44_46.tif 44_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 44_47.tif 44_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 44_48.tif 44_48_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 45_44.tif 45_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 45_45.tif 45_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 45_46.tif 45_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 45_47.tif 45_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 45_48.tif 45_48_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 46_44.tif 46_44_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 46_45.tif 46_45_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 46_46.tif 46_46_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 46_47.tif 46_47_noborder.tif
gdal_translate -a_nodata 111 -co tfw=yes 47_45.tif 47_45_noborder.tif

8.3. Metoda 3 (automată)

Folosiți un limbaj de programare pentru a genera sau afișa pe ecran textul fișierului .bat. Orice limbaj este bun. În continuare voi prezenta un exemplu de script scris în PHP. Fișierul rezultat seamănă cu cel produs folosind OpenOffice, doar ordinea comenzilor fiind puțin diferita (în loc să avem 40 de rînduri cu prima comandă, urmate de încă 40 cu a doua comandă, cele două comenzi vor alterna de 40 de ori la rind pentru fiecare foaie de hartă).

<?php

//tablou ce contine numele foilor de harta. Numele ar putea fi citite direct dintr-un folder, dintr-un fisier text sau o baza de date
$planse = array ("38_45", "38_46", "39_45", "39_46", "39_47", "40_44", "40_45", "40_46", "40_47", "40_48", "41_44", "41_45", "41_46", "41_47", "41_48", "42_44", "42_45", "42_46", "42_47", "42_48", "43_44", "43_45", "43_46", "43_47", "43_48", "44_44", "44_45", "44_46", "44_47", "44_48", "45_44", "45_45", "45_46", "45_47", "45_48", "46_44", "46_45", "46_46", "46_47", "47_45");

//variabila care va stoca continutul fisierului
$continut_fisier = "";

//parcurgerea intr-o bucla a tuturor numelor
foreach ($planse as $plansa)
{
	//crearea comenzilor de preprocesare a datelor
	$continut_fisier .= 'gdal_rasterize -b 1 -b 2 -b 3 -burn 111 -burn 111 -burn 111 -i -l '.$plansa.' '.$plansa.'.shp '.$plansa.'.tif'."\n";
	$continut_fisier .= 'gdal_translate -a_nodata 111 -co tfw=yes '.$plansa.' '.$plansa.'_noborder.tif'."\n";
}

//scriere fisier .bat pe disc
$f = fopen("procesare_harti.bat", "w");
fwrite($f, $continut_fisier);
fclose($f);

?>

8.4. Completare

Indiferent de metoda aleasă, completați fișierul .bat cu ultimele două linii:

gdalwarp -of Gtiff -co tfw=yes -co 38_45_noborder.tif 38_46_noborder.tif 39_45_noborder.tif 39_46_noborder.tif 39_47_noborder.tif 40_44_noborder.tif 40_45_noborder.tif 40_46_noborder.tif 40_47_noborder.tif 40_48_noborder.tif 41_44_noborder.tif 41_45_noborder.tif 41_46_noborder.tif 41_47_noborder.tif 41_48_noborder.tif 42_44_noborder.tif 42_45_noborder.tif 42_46_noborder.tif 42_47_noborder.tif 42_48_noborder.tif 43_44_noborder.tif 43_45_noborder.tif 43_46_noborder.tif 43_47_noborder.tif 43_48_noborder.tif 44_44_noborder.tif 44_45_noborder.tif 44_46_noborder.tif 44_47_noborder.tif 44_48_noborder.tif 45_44_noborder.tif 45_45_noborder.tif 45_46_noborder.tif 45_47_noborder.tif 45_48_noborder.tif 46_44_noborder.tif 46_45_noborder.tif 46_46_noborder.tif 46_47_noborder.tif 47_45_noborder.tif merge.tif
gdal_translate -of ECW -co "LARGE_OK=YES" merge.tif merge.ecw

Salvați modificările făcute și deschideți o instanță FWTools Shell. Navigați pînă în folderul unde ați descărcat foile de hartă, scrie-ți numele fișierului .bat în linia de comandă și apăsați tasta Enter. Seria de comenzi va fi rulată automat iar, în final, la capătul a aproximativ o oră de procesare, fișierul .ecw se va găsi stocat pe calculatorul propriu.

Concluzii

Pe parcursul acestui tutorial am încercat să vă familiarizez cu comenzile GDAL/OGR ce pot fi folosite pentru preprocesarea, mozaicarea și compresia datelor geospatiale raster. Cu foarte puțin efort, metodologia prezentată aici poate fi adaptată pentru un alt set de date. Pentru orice sugestii, completări sau nelămuriri vă așteptăm pe lista de discuții.