Transcript Interfacing Spatial Statistics sofware and GIS: a case study

Interfacing Spatial
Statistics sofware and
GIS: a case study
TerraLib: Open source GIS library

Data management


Functions


All of data (spatial + attributes) is in
database
Spatial statistics, Image Processing,
Map Algebra
Innovation

Based on state-of-the-art techniques
 Same timing as similar commercial
products

Web-based co-operative development

http://www.terralib.org
Operational Vision of TerraLib
DBMS
TerraLib
Geographic
Application
Spatial
Operations
API for
Spatial
Operations
Spatial
Operations
Access
Oracle
Spatial
MySQL
Postgre
SQL
TerraLib  MapObjects + ArcSDE + cell spaces + spatio-temporal models
TerraLib applications

Cadastral Mapping


Public Health


Indicators of social exclusion in innercity areas
Land-use change modelling


Spatial statistical tools for
epidemiology and health services
Social Exclusion


Improving urban management of large
Brazilian cities
Spatio-temporal models of
deforestation in Amazonia
Emergency action planning

Oil refineries and pipelines (Petrobras)
TerraLib Structure
Java Interface
COM Interface
OGIS Services
C++ Interface
Functions
kernel
Visualization
Controls
Spatio-Temporal
Data Structures
File and DBMS
Access
I/O Drivers
External
Files
DBMS
TerraView


TerraLib client
Connects to DBMS


Visualization with ‘views’ and ‘themes’



Images, Grids, vectors and tables
Brushing between graphics and maps
Analysis


Oracle, PostGIS, MySQL and SQLServer
Geocoding, ESDA, Clustering (more to come)
I/O


SHP, MIF, and SPRING vectors
TIFF and GEOTIFF rasters
Creating a MySQL spatial database in TerraView

File...Open Database menu
Importing GIS files

File...Import data



Select the rio_inf_mort.shp (Rio infant mortality data)
Use OBJET_ID_1 as the column
Use Projection... Lat/Long and Datum SAD69
Databases, layers, views and themes

A database has many
layers

This is a storage concept

Visualization is based on
views and themes

A view is a set of spatial
objects from a layer

A theme is a way to display
these objects
Changing the legend of a theme

Theme...Change legend

Change the legend to
display the NEO_RATE


Quantiles
Five slices
Resulting TerraView interface
Producing a scatterplot


Theme...Graphic Parameters
Selection a DISPERSION graphic


X – TOTBIRTH
Y – NEO_RATE
Brushing

Operation...TileWindows


Shows both graphics and map
Brush between the two and the table
R – TerraLib integration
Strong coupling mechanism
R – TerraLib integration

Enables R to access a TerraLib database

Results in R are incorporated in a GIS database directly

Data is analysed in R using packages such as gstat or
geoR

Results are visualized in TerraView

Integration is performed by the aRT API
aRT: API R-TerraLib

aRT is an R package that provides the integration
between the software R and the GIS classes TerraLib.

Developed by UFPR, Brazil

http://www.est.ufpr.br/aRT/
aRT: API R-TerraLib

Classes to manipulate data and functions in TerraLib

aRT


aRTdb


Creation and access to a database
aRTlayer


Enables connection to the DBMS for database administration
Enables manipulation of layers in aRT
aRTtheme

Enables visualização of themes in TerraView
Integration R - TerraLib
R graphics
TerraView map
(Neto et al., 2005)
R-TerraLib integration: example 1

Load the geoR and arT packages into R manually

geoR is a library for advanced geostatistics (also
developed by Paulo Ribeiro, UFPR)

Get the geoR and arT packages from the course
homepage

Expand the zip files and copy them to the directory
C:\Program Files\R\rw2011\library
R-TerraLib integration: example 1
##
## loading packages
##
require(geoR)
require(aRT)
##
## Exemple 1: data for parana state
##
data(parana)
points(parana, bor=borders)
A new type: geodata
> class (parana)
[1] "geodata"
> parana
$coords
east
north
[1,] 402.9529 164.52841
[2,] 501.7049 428.77100
[3,] 556.3262 445.27065
[4,] 573.4043 447.04177
[5,] 702.4228 272.29590
[6,] 668.5442 261.67070
[7,] 435.8477 286.54044
[8,] 434.0125 317.90506
[9,] 432.4622 288.37001
A new type: geodata
$data
[1]
[8]
[15]
[22]
[29]
[36]
[43]
[50]
[57]
[64]
[71]
306.09
282.07
237.87
351.73
199.91
212.50
321.69
302.13
344.59
236.25
165.46
200.88
319.12
222.87
277.92
167.00
242.08
208.89
335.41
366.10
228.82
320.31
167.07
244.67
263.10
323.08
182.88
247.79
238.62
330.87
201.84
258.12
232.80
162.77
233.31
236.91
253.32
197.10
187.27
248.76
329.49
218.27
232.17
266.27
163.57
224.46
247.01
315.33
257.50
222.54
193.48
262.81
200.38
248.17
301.10
178.61
206.12
240.58
379.94
205.16
313.60
240.48
365.88
229.40
240.66
244.78
301.54
248.99
304.28
197.09
224.07
269.92
265.56
359.08
235.07
184.59
248.57
A new type: geodata
$borders
[1,]
[2,]
[3,]
[4,]
[5,]
[6,]
[7,]
[8,]
[9,]
[10,]
[11,]
[12,]
[13,]
east
670.2049
663.7187
656.0667
649.9714
642.6346
635.2717
628.7421
622.7559
617.2708
610.4570
604.4217
600.2620
592.6210
north
111.7610
107.0510
105.2420
100.7915
97.9930
101.1545
105.5405
110.7495
116.3970
120.5795
119.6670
121.8980
119.9675
A new type: geodata
$loci.paper
[,1]
[1,] 300.336
[2,] 647.755
[3,] 361.764
[4,] 410.100
[,2]
484.453
317.463
438.781
260.373
R-TerraLib integration: example 1
##
## connection to MySQL
##
> con <- openConn(user=“gilberto")
> con
Object of class aRTconn
DBMS: "MySQL"
User: "gilberto"
Password: ""
Port: 3306
Host: "localhost"
Cleaning the database

If the database exists, clean it
> showDbs(con)
[1] "ifgi"
"intro"
"test"
"mysql"
"parana"
> if(any(showDbs(con)=="parana"))
deleteDb(con, "parana", force=T)
Deleting database 'parana' ... yes
Creating a new database
## creating a new database
> pr = createDb(con, "parana")
Connecting with database 'parana' ... no
Creating database 'parana' ... yes
Creating conceptual model of database
'parana' ... yes
Loading layer set of database 'parana' ...
yes
Loading view set of database 'parana' ... yes
Creating a layer in the database
##
## Creating a layer to hold the data
##
 l_data <- createLayer(pr, "data")
Building projection to layer 'data' ... yes
Creating layer 'data' ... yes
Loading points into the database
## loading points
> artcoords <- .gr2aRTpoints(parana)
# adding points to a layer
> addPoints(l_data, artcoords)
Converting points to TerraLib format ... yes
Adding 143 points to layer 'data' ... yes
Reloading tables of layer 'data' ... yes
# adding a table to the layer
> t_data = createTable(l_data, "t_data",gen=T)
Creating static table 't_data' on layer 'data'
... yes
Creating link ids ... yes
Plotting the data
> plot(l_data)
Inserting attributes on the GIS table
> artdata <- data.frame(.gr2aRTattributes(parana))
>
> artdata[1:5,]
id
data
1 1 306.09
2 2 200.88
3 3 167.07
4 4 162.77
5 5 163.57
>
> names(artdata)
[1] "id"
"data"
> createColumn(t_data, "data")
Checking for column 'data' in table 't_data' ... no
Creating column 'data' in table 't_data' ... yes
> updateColumns(t_data, artdata)
Converting 2 attributes to TerraLib format ... yes
Converting 143 rows to TerraLib format ... yes
Updating columns of table 't_dados' ... yes
Creating views and themes for TerraView
## criando vistas e temas automaticas para o TV
(opcional)
> th_data <- createTheme(l_data, “stations", view =
"view")
Checking for theme ‘stations' in layer 'parana' ...
no
Creating theme ‘stations' on layer 'dados' ... yes
Checking for view 'view' in database 'parana' ... no
Creating view 'view' ... yes
Inserting view 'view' in database 'parana' ... yes
Checking tables of theme ‘stations' ... yes
Saving theme ‘stations' ... yes
Building collection of theme ‘stations' ... yes
> setVisual(tema_dados, visualPoints(pch=22,
color="black", size=5))
Creating a layer to store borders
> artpols <- .gr2aRTpolygons(parana)
>
> l_pol<-createLayer(pr, “borders")
Building projection to layer ‘borders' ... yes
Creating layer ‘borders' ... yes
> addPolygons(l_pol, artpols)
Converting polygons to TerraLib format ... yes
Adding 1 polygons to layer ‘borders' ... yes
Reloading tables of layer ‘borders' ... yes
> createTable(l_pol, "t_pol")
Creating static table 't_pol' on layer ‘borders' ... yes
Creating link ids ... yes
Object of class aRTtable
Table: "t_pol"
Type: static
Layer: “borders"
Rows: 1
Attributes:
id: string[16] (key)
> plot(l_pol)
Creating a view and theme for the polygons
> tema_pol<-createTheme(l_pol, “borders",
view="view")
Checking for theme ‘borders' in layer 'parana' ...
no
Creating theme ‘borders' on layer ‘borders' ... yes
Checking for view 'view' in database 'parana' ...
yes
Checking tables of theme ‘borders' ... yes
Saving theme ‘borders' ... yes
Building collection of theme ‘borders' ... yes
> setVisual(tema_pol, visualPolygons())
Creating a grid for the interpolated surface
>
>
>
>
>
>
>
gx <- seq(50,800, by=10)
gy <- seq(-100,650, by=10)
loc0 <- expand.grid(gx,gy)
points(parana, bor=borders)
points(loc0, pch=".", col=2)
loc1 <- polygrid(loc0, bor=parana$bor)
points(loc1, pch="+", col=4)
Kriging using geoS: model fitting
> # Kriging using geoS
> # Step 1 – fitting the variogram
> ml <- likfit(parana, trend="1st", ini=c(1000,
100))
--------------------------------------------------likfit: likelihood maximisation using the function
optim.
likfit: WARNING: This step can be time demanding!
--------------------------------------------------likfit: end of numerical maximisation.
Kriging using geoS: interpolation
> kc <- krige.control(trend.d="1st", trend.l="1st",
obj=ml)
> kc <- krige.conv(parana, loc=loc0, krige= KC,
bor=parana$borders)
krige.conv: results will be returned only for
prediction locations inside the borders
krige.conv: model with mean given by a 1st order
polynomial on the coordinates
krige.conv: Kriging performed using global
neighbourhood
> save.image()
>
> image(kc, col=terrain.colors(15))
Exporting the grid to TerraLib
> # preparing object for aRT
> georpred <- .prepare.graph.kriging(locations=loc0,
borders=parana$borders, values=kc$pred)
> artpred <- .gr2aRTraster(georpred)
> # creating a layer in TerraLib
> l_pred <- new("aRTlayer", pr, layer="pred",
create=T)
Building projection to layer 'pred' ... yes
Creating layer 'pred' ... yes
> addRaster(l_pred, artpred)
Adding raster data to layer 'pred' ... yes
Reloading tables of layer 'pred' ... yes
Creating a theme and a view for the grid
> th<-createTheme(l_pred, "raster",
view="view")
Checking for theme 'raster' in layer
'parana' ... no
Creating theme 'raster' on layer 'pred' ...
yes
Checking for view 'view' in database
'parana' ... yes
Checking tables of theme 'raster' ... yes
Saving theme 'raster' ... yes
>
> setVisual(th, visualRaster(), mode="r")
Plotting the data in R
>
>
>
>
# plotting the data in R
plot(l_pred, col=terrain.colors(15))
plot(l_dados, add=T)
plot(l_pol, add=T)
Now, let’s see the data in TerraView


File...OpenDatabase...
Parana database is there!
Create Views and Themes (if needed)

Create View parana

Create themes




Borders
Data
Pred
Play with the data!
R data in TerraView