Bulk Interpolation using R Environment - PowerPoint Presentation

Slide1

Bulk Interpolation using R Environment

Jiří Kadlec – Aalto University, Finland Pavel Treml – Masaryk Water Research Institute and Charles University, Czech Republic20.9.2013 FOSS4G Conference - NottinghamSlide2

Field

s with Observations in time and spaceWhat, Where, When (irregular sampling)

Space, S

Time, T

Variables, V

s

t

V

i

D

c

“Where”

“What”

“When”

A data value

4.2

Praha

2013-09-20

Air Temperature (C)

Image created by CUAHSI, 2010Slide3

Type of data in this tutorial

Tens of variablesHundreds of stationsHundreds of observations per stationIncomplete dataObservations

station

longitudelatitudevariabletime

valuePraha14.43050.0156temperature

2013-03-31-0.2Praha14.43050.0156

temperature

2013-04-01

1.1

Praha

14.430

50.0156snow2013-04-01 3.0Brno16.25449.147temperature

2013-03-31-0.1Brno16.25449.147snow2013-03-0112.0

TASKFor each time-step, Examine how one orMore variables changeIn spaceSlide4

Interpolation

Deterministic methods

Geostatistical

methodsRepeated over many time stepsHow to automate?

?

time

value

lat

lon

?

Interpolation - Time

Interpolation - SpaceSlide5

the R Environment

Free statistical softwareWindows, Mac, LinuxCreate High-quality graphicsMany input data formatsText file Vector data (shapefile)Raster data (grid)Many output data formatsPictureText file

Raster, vectorScripting language for automating repeated tasksSlide6

Case study: Maps of air

temperature andsnow, Czech RepublicYear 2013 (365 maps)Require identical color-scaleNeed to show rivers, boundariesNeed to show point values (to assess interpolation)

Load Data sets

Select next time Select matching observations

Run Interpolation

Create mapcartography

More times?

END

START

NO

YESSlide7

Load Data Sets – Raster (SRTM elevation)

Library(“sp”)library(“raster")

library(“rcolorBrewer")srtm <- raster("srtm_utm.asc")colors=

brewer.pal(6, "YlOrRd")intervals = c(0, 250, 500, 750, 1000, 1600)spplot(

srtm_proj, col.regions=colors, at=intervals)

Select packages

With needed functions

Color ramp setting

Raster:

Load from local file

Visualize raster using

spplot

Note:

You can get DEM for your area in R using:Slide8

Load Data Sets – Vector (boundaries, rivers)

library(“maptools")library(“sp")border = readShapeLines(“border.shp")border_layer = list("

sp.lines", border, lwd=2.0,col=“black")rivers = readShapeLines(“rivers.shp")

proj4string(rivers) <- CRS("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333")rivers_proj <- spTransform(rivers, CRS("+proj

=utm +zone=33"))river_layer = list("sp.lines", rivers_proj, lwd=1,col=“blue")

layout = list(river_layer, border_layer) spplot(srtm, sp.layout

=layout)

Select packages

With needed functions

Rivers:

Load

shapefile

,

Reproject vector fromKrovak to UTM system

Borders:

Load shapefile

Visualize

vector datasets

on top of rasterSlide9

Load Text File – Point Data (Stations)

At first, we can use a subset (for first date/time in the dataset)data = read.table(“data.txt”,

header = TRUE, sep = "\t", dec = ".")st = subset(data, DateTimeUTC == ‘2012-08-12’ &

VariableCode==“SNOW")Coordinates(st) = ~Long + Latproj4string(st) = CRS("+init

=epsg:4326")stations <- spTransform(st, CRS("+proj=utm +zone=33"))

stations_layer = list("sp.points", stations, pch=19, cex=1.0, col="black")labels = list("panel.text

",

coordinates(stations)[,

1],

coordinates(stations)[,

2],

labels=stations$value, col="black", font=1, pos=2)layout = list(stations_layer, labels)

Reproject from Lat/Lon to UTM systemSlide10

All Layers together (raster, vector, points)

(We use different color ramp in this example)Slide11

Interpolation: IDW and

Kriging MethodsSlide12

Interpolation: Covariate (Elevation)

model = lm(value ~ elev, data)intercept = model$coefficients[1]slope = model$coefficients[2]plot(value~elev, data)

abline(model)TEMP = a * ELEVATION + b

In our case temperature is often negativelycorrelated with elevationSlide13

Interpolation: Elevation as covariate

TEMP = a * ELEVATION + b

residuals

Instead of interpolating temperature directly,

we create grid using regression equation andwe only interpolate the residuals

Using TEMP = a * ELEVATION + b

Interpolated residuals

Combination (regression +

interpol

.

residualsSlide14

Color Breaks, Color Ramps

#user-defined color breakscolors = rev(brewer.pal(8, "YlGnBu"))brk <- c(-40,-35,-30,-25,-20,-15,-10,-5,0)plot(grid,breaks=brk,col=palette(colors)

One

color ramp

and set of

user-defined color breaks easily re-used for different grids

RColorBrewer

packages provides

Pre-defined color rampsSlide15

Example Final Map: One time step

Combines the previous steps …Snow depth (cm): 2013-01-01Slide16

Saving map to file

figure = “2012-02-07.jpg”png(filename = figure, width = 1500, height = 1000, pointsize = 25, quality = 100, bg = "white", res = 150)DO THE PLOT COMMANDS HEREdev.off()

Set image sizeSet image resolution

Other options (margins, multiple maps in one image)

PNG

picture

JPEG

picture

PDF

file

WMF

file

.asc ascii grid file (for GIS softwares)……..Slide17

Bulk Interpolation: “FOR” loop

(multiple time steps)for (j in 1:length(timesteps)) { # select subset of observations # run interpolation # create

map # export map to file (picture, pdf, …)}

Schema of the LoopLoad Data sets

Select next time

Select matching observations

Run Interpolation

Create map

cartography

More times?

END

START

NO

YESSlide18

Final Map: Multiple time stepsSlide19

R as Compared with Desktop GIS

R Statistical Software Environment

Desktop GISMaps are highly interactiveCreate small number of maps with

graphical user interfaceAutomated cartography (label placement…)Map-Layout, model builder toolsArcGIS, QGIS, SAGA, MapWindow

,…IDV, Panoply, … (-) Maps are static

(-) Some commands counter-intuitive for novice user(+) Create very large number of mapswith same cartographic symbology

(+)

Task

is easy to

automate

and

reproducible(+) Good documentation and large user baseDesktop GIS: project file (.mxd, .qpj, …)R: script (.R)Slide20

References

Bivand, Roger S., Edzer J. Pebesma, and Virgilio Gómez Rubio. Applied spatial data: analysis with R. Springer, 2008.Hengl, Tomislav. A practical guide to geostatistical mapping of environmental variables. Vol. 140. No. 4. 2007.

BOOKS

www.spatial-analyst.net

www.r-project.org

http://

hydrodata.info/api/ (

source of data in our demos)

WEBSITES

R-PACKAGES USED

RColorBrewer

,

maptools

,

rgdal

,

sp

, rasterSlide21

Try the scripts yourself!

hydrodata.info/R/Sample scripts for bulk interpolationThis presentationSource data { central Europe meteorological observations }