Interfacing Spatial Statistics sofware and GIS: a case study

1
Interfacing Spatial Statistics sofware and GIS a
case study
2
TerraLib Open source GIS library

• Data management
• All of data (spatial attributes) is in database
• Functions
• 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

3
Operational Vision of TerraLib
TerraLib ? MapObjects ArcSDE cell spaces
spatio-temporal models
4
TerraLib applications

• Cadastral Mapping
• Improving urban management of large Brazilian
  cities
• Public Health
• Spatial statistical tools for epidemiology and
  health services
• Social Exclusion
• Indicators of social exclusion in inner-city
  areas
• Land-use change modelling
• Spatio-temporal models of deforestation in
  Amazonia
• Emergency action planning
• Oil refineries and pipelines (Petrobras)

5
TerraLib Structure
Java Interface
OGIS Services
COM Interface
C Interface
Functions
kernel
Spatio-Temporal Data Structures
File and DBMS Access
Visualization Controls
I/O Drivers
DBMS
External Files
6
TerraView

• TerraLib client
• Connects to DBMS
• Oracle, PostGIS, MySQL and SQLServer
• Visualization with views and themes
• Images, Grids, vectors and tables
• Brushing between graphics and maps
• Analysis
• Geocoding, ESDA, Clustering (more to come)
• I/O
• SHP, MIF, and SPRING vectors
• TIFF and GEOTIFF rasters

7
Creating a MySQL spatial database in TerraView

• File...Open Database menu

8
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

9
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

10
Changing the legend of a theme

• Theme...Change legend
• Change the legend to display the NEO_RATE
• Quantiles
• Five slices

11
Resulting TerraView interface
12
Producing a scatterplot

• Theme...Graphic Parameters
• Selection a DISPERSION graphic
• X TOTBIRTH
• Y NEO_RATE

13
Brushing

• Operation...TileWindows
• Shows both graphics and map
• Brush between the two and the table

14
R TerraLib integration

• Strong coupling mechanism

15
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

16
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/

17
aRT API R-TerraLib

• Classes to manipulate data and functions in
  TerraLib
• aRT
• Enables connection to the DBMS for database
  administration
• aRTdb
• Creation and access to a database
• aRTlayer
• Enables manipulation of layers in aRT
• aRTtheme
• Enables visualização of themes in TerraView

18
Integration R - TerraLib
19
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

20
R-TerraLib integration example 1

• loading packages
• require(geoR)
• require(aRT)
• Exemple 1 data for parana state
• data(parana)
• points(parana, borborders)

21
A new type geodata

• gt class (parana)
• 1 "geodata"
• gt 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

22
A new type geodata

• data
• 1 306.09 200.88 167.07 162.77 163.57 178.61
  301.54
• 8 282.07 319.12 244.67 233.31 224.46 206.12
  248.99
• 15 237.87 222.87 263.10 236.91 247.01 240.58
  304.28
• 22 351.73 277.92 323.08 253.32 315.33 379.94
  197.09
• 29 199.91 167.00 182.88 197.10 257.50 205.16
  224.07
• 36 212.50 242.08 247.79 187.27 222.54 313.60
  269.92
• 43 321.69 208.89 238.62 248.76 193.48 240.48
  265.56
• 50 302.13 335.41 330.87 329.49 262.81 365.88
  359.08
• 57 344.59 366.10 201.84 218.27 200.38 229.40
  235.07
• 64 236.25 228.82 258.12 232.17 248.17 240.66
  184.59
• 71 165.46 320.31 232.80 266.27 301.10 244.78
  248.57

23
A new type geodata

• borders
• east north
• 1, 670.2049 111.7610
• 2, 663.7187 107.0510
• 3, 656.0667 105.2420
• 4, 649.9714 100.7915
• 5, 642.6346 97.9930
• 6, 635.2717 101.1545
• 7, 628.7421 105.5405
• 8, 622.7559 110.7495
• 9, 617.2708 116.3970
• 10, 610.4570 120.5795
• 11, 604.4217 119.6670
• 12, 600.2620 121.8980
• 13, 592.6210 119.9675

24
A new type geodata

• loci.paper
• ,1 ,2
• 1, 300.336 484.453
• 2, 647.755 317.463
• 3, 361.764 438.781
• 4, 410.100 260.373

25
R-TerraLib integration example 1

• connection to MySQL
• gt con lt- openConn(usergilberto")
• gt con
• Object of class aRTconn
• DBMS "MySQL"
• User "gilberto"
• Password ""
• Port 3306
• Host "localhost"

26
Cleaning the database

• If the database exists, clean it
• gt showDbs(con)
• 1 "ifgi" "intro" "mysql" "parana" "test"
• gt if(any(showDbs(con)"parana")) deleteDb(con,
  "parana", forceT)
• Deleting database 'parana' ... yes

27
Creating a new database

• creating a new database
• gt 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

28
Creating a layer in the database

• Creating a layer to hold the data
• l_data lt- createLayer(pr, "data")
• Building projection to layer 'data' ... yes
• Creating layer 'data' ... yes

29
Loading points into the database

• loading points
• gt artcoords lt- .gr2aRTpoints(parana)
• adding points to a layer
• gt 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
• gt t_data createTable(l_data, "t_data",genT)
• Creating static table 't_data' on layer 'data'
  ... yes
• Creating link ids ... yes

30
Plotting the data

• gt plot(l_data)

31
Inserting attributes on the GIS table

• gt artdata lt- data.frame(.gr2aRTattributes(parana))
  
• gt
• gt artdata15,
• id data
• 1 1 306.09
• 2 2 200.88
• 3 3 167.07
• 4 4 162.77
• 5 5 163.57
• gt
• gt names(artdata)
• 1 "id" "data"
• gt createColumn(t_data, "data")
• Checking for column 'data' in table 't_data' ...
  no
• Creating column 'data' in table 't_data' ... yes
• gt 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

32
Creating views and themes for TerraView

• criando vistas e temas automaticas para o TV
  (opcional)
• gt th_data lt- 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
• gt setVisual(tema_dados, visualPoints(pch22,
  color"black", size5))

33
Creating a layer to store borders

• gt artpols lt- .gr2aRTpolygons(parana)
• gt
• gt l_pollt-createLayer(pr, borders")
• Building projection to layer borders' ... yes
• Creating layer borders' ... yes
• gt addPolygons(l_pol, artpols)
• Converting polygons to TerraLib format ... yes
• Adding 1 polygons to layer borders' ... yes
• Reloading tables of layer borders' ... yes
• gt 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 string16 (key)

34
Creating a view and theme for the polygons

• gt tema_pollt-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
• gt setVisual(tema_pol, visualPolygons())

35
Creating a grid for the interpolated surface

• gt gx lt- seq(50,800, by10)
• gt gy lt- seq(-100,650, by10)
• gt loc0 lt- expand.grid(gx,gy)
• gt points(parana, borborders)
• gt points(loc0, pch".", col2)
• gt loc1 lt- polygrid(loc0, borparanabor)
• gt points(loc1, pch"", col4)

36
Kriging using geoS model fitting

• gt Kriging using geoS
• gt Step 1 fitting the variogram
• gt ml lt- likfit(parana, trend"1st", inic(1000,
  100))
• --------------------------------------------------
  -
• likfit likelihood maximisation using the
  function optim.
• likfit WARNING This step can be time demanding!
• --------------------------------------------------
  -
• likfit end of numerical maximisation.

37
Kriging using geoS interpolation

• gt kc lt- krige.control(trend.d"1st",
  trend.l"1st", objml)
• gt kc lt- krige.conv(parana, locloc0, krige KC,
  borparanaborders)
• 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
• gt save.image()
• gt
• gt image(kc, colterrain.colors(15))

38
Exporting the grid to TerraLib

• gt preparing object for aRT
• gt georpred lt- .prepare.graph.kriging(locationsloc
  0, bordersparanaborders, valueskcpred)
• gt artpred lt- .gr2aRTraster(georpred)
• gt creating a layer in TerraLib
• gt l_pred lt- new("aRTlayer", pr, layer"pred",
  createT)
• Building projection to layer 'pred' ... yes
• Creating layer 'pred' ... yes
• gt addRaster(l_pred, artpred)
• Adding raster data to layer 'pred' ... yes
• Reloading tables of layer 'pred' ... yes

39
Creating a theme and a view for the grid

• gt thlt-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
• gt
• gt setVisual(th, visualRaster(), mode"r")

40
Plotting the data in R

• gt plotting the data in R
• gt plot(l_pred, colterrain.colors(15))
• gt plot(l_dados, addT)
• gt plot(l_pol, addT)

41
Now, lets see the data in TerraView

• File...OpenDatabase...
• Parana database is there!

42
Create Views and Themes (if needed)

• Create View parana
• Create themes
• Borders
• Data
• Pred
• Play with the data!

43
R data in TerraView