Harta repartiției teritoriale a cantităților medii de precipitații: Residual Kriging (metoda simplă). Aplicație ILWIS

25.10.2014 Tutorialul de fața folosește o versiunea mai veche a softului prezentat. Va rugăm consultați versiunea actualizată sau contactați autorul.

Distribuția spațială a cantităților medii de precipitații este puternic influențată, la fel ca în cazul valorilor de temperatură , de altitudine. Așadar, în realizarea hărții cu repartiția teritorială valorilor medii de precipitații, trebuie întotdeauna testată intensitatea legăturii dintre cele două variabile, iar atunci când este cazul altimetria trebuie luată în considerare. Pentru realizarea hărții de precipitații se va folosi programul ILWIS Open 3.7 utilizând Residual Kriging (o metodă de interpolare hibridă care poate include în procesul de spațializare una sau mai multe variabile auxiliare).

Etapele realizării harții sunt:

1. Aflarea relației dintre parametrul analizat și variabila explicativă (altitudine) prin aplicarea regresiei simple liniare;
2. Calcularea reziduurilor regresiei (diferența dintre valorile estimate prin regresie și cele măsurate) pentru fiecare locație;
3. Modelarea spațială a valorilor de precipitații prin aplicarea regresiei spațiale utilizând modelul numeric al terenului¹;
4. Interpolarea reziduurilor folosind metoda Ordinary Kriging;
5. Însumarea suprafeței modelate (punctul 3) cu cea a reziduurilor (punctul4).

Prin aplicarea ecuației de regresie dintre varibilă independentă (altitudinea) și variabila dependentă (valorile de precipitații) rezultă un coeficient de determinare de R² =0,79, mai exact 79 % din variabilitatea spațială a valorilor medii de precipitații este explicată de altimetrie.

Pentru calcularea coeficientului de corelație datele au fost aduse in formatul ILWIS:

		import shape(..\grided_precipitation\r\pp_61_90.shp, pp_61_90)	
	

iar altitudinea fiecărei stații meteorologice a fost extrasă (în tabelul pp_61_90_corelatie) din modelul numeric al terenului utilizându-se comanda:

		dem=MAPVALUE(dem,Coordinate)	
	

Coeficientului de corelație (r=0.891) se calculează utilizând commanda:

		corelatie=corr(dem,pp_august)	
	

iar prin ridicarea la pătrat a acestuia se obține coeficientul de determinare (R² = 0.794):

		r2=(corelatie^2)	
	

Coeficientul de corelație (r) reprezintă un instrument statistic de măsură a gradului de dependență dintre doua variabile (precipitații și altitudine în cazul nostru). Coeficientul de corelație variază intre -1 si +1. Cu cât coeficientul de corelație este mai apropiat de valoarea 1, cu atât relația dintre cele două variabile analizate este mai puternică. La valori negative, coeficientul indica o corelație inversă, iar atunci când valoarea acestuie este apropiată de 0, indica lipsa unei corelații.

Coeficientul de determinare (R²) este dat de pătratul coeficientului de corelație (r). Acesta indică cât de bine dreapta de regresie aproximează norul de puncte. R² reprezintă măsura în care variabila independentă, X (altitudinea), explică variația variabilei dependente Y (precipitațiile). O valoare dezirabilă a coeficientului de determinare este una cat mai apropiată de 1.

Pentru aflarea dreptei de regresie care aproximează tendința norului de puncte între cele două variabile s-a realizat graficul dintre valorile precipitatțiilor și altidudinilor stațiilor (Figura 1). Această relație este dată de ecuația y = 48.117+0.043*x.

Figura 1. Dreapta de regresie.

Calculare reziduurilor (erorilor regresiei) se obține prin diferența (residuu=pp_august-pp_august_regress) dinte valorile măsurate și cele estimate obținute utilizând relația (Figura 2):

		pp_august_regress=48.117+0.043*dem	
	

Figura 2. Tabelul pp_61_90_corelatie.

Gridul cu modelarea spațială a valorilor de precipitații este obținut prin aplicarea regresiei spațiale utilizând modelul numeric al terenului (dem) utilizând comanda:

		pp_august_regress=48.117+0.043*dem
	

Suprafața erorilor regresiei rezultă prin interpolarea câmpului residuu din tabelul pp_61_90_corelatie utilizând metoda Ordinary Kriging (krigingul normal/ordinar) care foloseste combinațiile liniare ponderate ale datelor (valorile precipitațiilor înregistrate la stațiile meterologice) utilizate în interpolare (Johnston et al., 2001):

Valoarea ponderii fiecărui punct folosit în interpolare rezultă în urma rezolvării unui sistem de ecuații utilizând calculul matricial de tipul (Hendrikse,2000):

Matricea din partea stângă depinde de distanța dintre punctele cu mărimi cunoscute și modelul semivariogramei ales. Matricea din partea dreaptă este dată de distanțele dintre locația aleasă pentru estimare și punctele cu măsurători. Aceasta va fi calculată pentru fiecare locație pentru care se realizează estimarea, în cazul nostru pentru fiecare pixel care formează gridul final.

Pentru calcularea parametrilor (nugget, sill, range), care sunt utilizați în rezolvarea sistemului de ecuații amintit mai sus, este necesar a se alege un model teoretic de semivariogramă potrivit datelor ce vor fi interpolate (reziduurile în cazul nostru). Din păcate ILWIS nu oferă un instrumental automat (intuitiv) în alegerea acestui model, pentru realizarea acestei etape utilizându-se programul Vesper 1.6 (Minasny, et al., 2002). Pentru aceste date s-a ales un model Exponențial cu parametrii nugget=9,15, sill=83,83 și range=62643 (Figura 3).

Figura 3. Modelul variogramei (Vesper 1.6).

Harta cu repartiția teritorială a reziduurilor regresiei (Figura 4) este obținută utilizând comanda:

		pp_residuu.mpr{dom=value;vr=-84.942:80.638:0.000} = MapKrigingOrdinary(pp_residuu.residuu,dem.grf,Exponential(9.15,83.83,62643),600000.000000,plane,0,1,60,no)
	

Figura 4. Harta reziduurilor.

Prin însumarea gridului reziduurilor cu cel al estimarii prin regresie se obține harta finală cu repartiția teritorială a cantitatților medii de precipitații (1961-1990) pentru luna august (figura 5):

		pp_august=pp_august_regress+pp_residuu
	

Figura 5. Repartiția teritorială a cantităților medii de precipitații, august 1961-1990.

Pentru validarea modelului de spațializare se comparară valorilor măsurate cu cele estimate calculându-se coeficientul de corelație r  care în cazul de față are o valoare de 0,995, valoarea 1 semnificând o perfectă corelație între cele două seturi de date analizate. Pentru a se extrage valorile estimate din gridul precipitațiilor medii lunare se folosește comanda:

		grid_estimat = mapvalue(pp_august, Coordinate)
	

iar pentru compararea celor două seturi de date (estimate vs. măsurat):

		pearson = CORR(pp_august, grid_estimat)
	

Bibliografie:

J. Hendrikse. Geostatistics in ILWIS. INTERNATIONAL ARCHIVES OF PHOTOGRAMMETRY AND REMOTE SENSING, 33(B4/1; PART 4):365-375, 2000.

K. Johnston, J.M. Ver Hoef, K. Krivoruchko, and N. Lucas. Using ArcGIS geostatistical analyst. Esri New York, 2001.

B. Minasny, AB McBratney, and BM Whelan. Vesper (variogram estimation and spatial prediction plus error) version 1.6. Aust. Ctr. for Precision Agric., Sydney, 2002.

¹ SRTM  1 km ² (http://www2.jpl.nasa.gov/srtm/)