Use of SAGA GIS for spatial interpolation (kriging) - DMCSEE

1 Use of saga GIS for spatial interpolation ( kriging ) Technical instructions Prepared for 1st DMCSEE -TCP training Budapest, 2-5 February 2010 Environmental Agency of Slovenia DELIVERABLE SUMMARY PROJECT INFORMATION Project acronym: DMCSEE Project title: Drought Management Centre for South East Europe Contract number: SEE/A/091 Starting date: 1. 4. 2009 Ending date: 31. 3. 2012 Project WEB site address: Lead partner organisation: Environmental Agency of the Republic of Slovenia Name of representative: dr. Silvo lebir, director Project manager: dr. GregorGregori E-mail: Telephone number: +386 (0)1 478 40 65 DELIVERABLE INFORMATION Title of the deliverable: saga kriging training manual WP/activity related to the deliverable: WP5/ Type (internal or restricted or public): Public Location (if relevant): WP leader: OMSZ Activity leader: Participating partner(s): OMSZ, EARS Author: Compiled by EARS E-mail: Contact: Telephone number: DELIVERY DEADLINES Contractual date of delivery to the JTS: Actual date of delivery to JTS: INTRODUCTION saga (System for Automated Geoscientific Analyses) is free open-source software, designed for implementation of spatial algorithms and geospatial visualization.

2 It was designed by University of Goettingen and further developed by University of Hamburg (see ). Although kriging algorithms, implemented within saga , lack some of features (flexible variogram fitting, cross validation) it can be used as tool for spatial interpolation of meteorological data, at least as first approximation. This document describes in practical manner steps to be taken in order to obtain maps of meteorological variables from point values ( station measurements). As example, interpolation of Standardized Precipitation Index (SPI) is presented for Slovenia, however, procedure is not limited to this variable and can be implemented for any spatial unit and any cartesian coordinate system. 1. PREPARATION OF DATA saga uses Shapefile format for vector data (points, lines, polygons). Therefore, it would be easiest if measurement data is prepared as shapefile that includes measurement values as point attributes. Since many database applications don t support data output in shapefile format, easiest approach would be to prepare data in plain ASCII file.

3 ASCII file should be organized in columns; among the columns, there have to be columns with X and Y coordinates; default is X as first column and Y as second column. Another requirement is for columns to be separated with single TAB character. If that is not the case (for instance, if your file is comma delimited), you can use simple ASCII editor (such as Notepad) to replace commas with TAB character (easiest way to do so is to store TAB character on Clipboard and paste it using Ctrl+V into Replace form in Notepad). First line should contain attribute names (pleas use X and Y for both coordinates) 2 Once your file is prepared (use extension .xyz that is default or .txt), import point data into saga : Modules Import/Export Shapes - Import Shapes From XYZ Enter column numbers for X and Y coordinates and file name. If import is successful, there should be Point data listed under Data tag carrying your file name. Recommendation: right-click on your data and save it as Shapefile.

4 Double-click on your data draws a map of all points: Map of all selected stations. Preparation of Variogram Modules Geostatistics points Semi-Variogram When window Semi-Variogram opens you put under: - Points : shape file with your point data - Attribute: you choose which attribute (in our case, which SPI length) will be used for variance calculations - Distance Increment: use a value that makes sense. In our case coordinates are in meter units; use at least 1000 (1 km) as distance increment. 3 Export to excel Save the table in saga : Under Data tag a table with name YourFileName[your_attribute]:Semi-varian ces should appear. Since saga doesn't provide mechanism to fit most frequently used variogram types, it is proposed to export calculated variances to MS Excel. Unfortunately it is not possible to export data via Clipboard (using Copy & Paste). You have to save variances as ASCII file.

5 Go under Data and open Semi-Variances table. With right mouse click you choose Save table as. Choose name and save variances in a .txt file. 4 Import to excel In Excel, go under Data Import external data Import data or simply choose Open file and choose Text files. You should import the Semi-Variances table that you have saved before in saga . 5 2. FITTING VARIOGRAM Variogram is a function describing dependence of variance with respect to distance. In order to apply kriging , variogram should be fitted to a simplified model. Usually there are three parameters that describe variogram to sufficient detail: nugget, sill and range. In Excel we will use formulas for spherical and exponential variogram model to fit variances to the data. Insert in the table formula for spherical variogram model as shown at the picture below: For example: if you put value of Range to cell M1, sill value to cell M2 and Nugget value to cell M3, the formula reads: =IF(A2<=$M$1;($M$2-$M$3)*(1,5*A2/$M$1 - 0,5*(A2/$M$1)^3)+$M$3;$M$2) where A2 is cell containing first distance class (imported from saga ).

6 IF(distance<=range;(sill-nugget)*(1,5*di stance/range - 0,5*(distance/range)^3)+nugget;sill) 6 In some cases exponential variogram model fit better: In the above specific example, Excel formula reads: =($M$2-$M$3)*(1 - EXP(-3*A2/$M$1)) +$M$3 Now you have formula only in first line. Select cell containing formula and drag formula over all distance classes. Make a graph: 1. series: X = Distance / Y = Variance 2. series: X = Distance / Y = Spherical or exponential model With changing Range, Sill and Nugget you can adjust model manually to fit the data. (sill-nugget)*(1 - EXP(-3*distance/range)) +nugget 7 3. kriging Ordinary kriging Modules Geostatistics- kriging Ordinary kriging (Global) First you have to decide to what raster you want interpolation to be performed. There are two options: - you can load an existing raster in same coordinate system as point coordinates have in the shapefile. In that case, spatial extent and size of pixel will be defined by that raster.

7 In that case, you should choose as target grid option grid system. - you can specify resolution (grid size) and coordinate extent (minimum and maximum values for X and Y coordinate) manually. In that case, you should choose as target grid option user defined. When window Ordinary kriging opens you put under: - Points: shape file - Attribute: you choose which SPI you want to be shown on the map - Options: untick Create variance grid for Target grid you choose grid system or user defined Variogram model Spherical model or exponential model Under Nugget, Sill and Range you put in numbers that you previously set in excel that fit the data. 8 By pressing "Okay" (and optionally choosing existing grid system) kriging process begins. It can take few minutes; in case of success, new grid appears under Data tag. You can plot a map by double click . Result: Map for SPI3 ( ) Note: resulting map has rectangular shape. For better orientation, further processing is needed (underlying with geographic maps, multiplying with mask; for latter, you can use modules > Grid calculus > grid calculator and apply a formula).

8 Universal kriging In some cases there is considerable correlation between analyzed variable and certain other spatial variable - such as altitude, east-west coordinate or others. In that case use of universal kriging is preferred: Is such correlation is expected then you have to prepare correlated variable as attribute in point measurements file. Show scatterplot: First calculate variances (as for ordinary kriging ).- Go under Data and open Semi-Variances table. With right mouse click you choose Show scatterplot 9 It can be seen that in this case X coordinate has significant correlation with SPI3 value: Therefore we will use universal kriging : Before that we have to prepare a grid with correlated variable; if correlated variable is altitude, you will need to import digital terrain layer in desired resolution. Modules Geostatistics- kriging Universal kriging (Global) When window Universal kriging opens you have to specify following information: - Points: shape file - Attribute: you choose which SPI you want to be shown on the map - Grids: enter pre-prepared grid layer containing only X coordinates - Options: untick Create variance grid for Target grid choose grid system Variogram model Spherical model Under Nugget, Sill and Range put in numbers that you previously set in Excel 10 Maps, resulting from universal kriging procedure, usually have visual impact of correlated variable.

9 In our case, impact of east-west gradient (which is obvious from general spatial pattern) is stronger when universal kriging is applied. Ordinary kriging ( ) Universal kriging ( ) 11 4. SMOOTHING OF MAPS In many cases when applying spatial interpolation , undesirable "bull eye" features appear. They result from isolated extreme values and are consequence of moderate spatial density of point measurements. Such a situation is shown in the picture below (difference between two SPI values of neighboring stations is large: -0,94 and 0,5). There are some ways how to reduce these features. Location of problematic spot The difference between two neighboring values is too large. 12 050100150200250300350 Semi Cerov ecMetlika rnomeljNov o mes toDvorKo evjes tationstotal rainfallJulyA ugus tSeptember1. Check if all the data for problematic stations is correct (no missing values in the data ..) 2. Erase station Look for stations nearby that are not included in shapefile and check their data.

10 For SPI calculations, criteria to include value is rather strict (long time series, nor interruptions or relocations). In case of isolated maxima, we can check values from stations not fulfilling this criteria - just as orientation. In the table below monthly sum of precipitation for some stations in the vicinity of "bull eye" are displayed (only stations Semi and Cerovec have appropriate data for SPI calculations and are originally included in the map). Data from other stations just help us determine that Semi really stands out so it won t be included in the map. station SPI3 July August September Semi -0,94 21,8 140,4 98,7 Cerovec 0,5 27,7 230,1 113 Metlika 28,5 220,4 98,4 rnomelj 25,6 266,1 122,3 Novo mesto 55,3 262,5 67,4 Dvor 36,4 241,6 123,9 Ko evje 50,1 323,9 136,8 3. Add data Instead of removing, we can calculate SPI for additional stations (even if they are not fulfilling the criteria) and add them to map in order to remove "bull eye" features.