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 http://earth.unibuc.ro/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.