In hazaeljones/geozoning: Zoning Methods for Spatial Data

  library(geozoning)
  library(sp)
  library(fields)

This vignette illustrates the zoning with corrections procedure on simulated data.

Generate a map

A map object is simulated with an exponential field and a variogram model. 450 points (default) are randomly allocated on a square field of size 1. Then 1936 points are kriged on a regular grid using inverse distance weighted interpolation. A Delaunay tesselation yields point neighborhood in the sense of Voronoi. For this purpose, we use the genMap function which is a wrapper of the randKmap function (see the documentation).

  seed=80
  map=genMap(DataObj=NULL,seed=seed,disp=FALSE,krig=2,Vmean=15,typeMod="Exp")

Display 2D map with three different views: first one=kriged data, second one=contour lines, third one=raw data.

  plotMap(map)

The map object contains all components required for zoning generation, evaluation and visualization. The data-related components include the kriged data, grid resolution, size in each dimension and map boundary. Kriged data are available as a SpatialPointsDataFrame object, as well as a matrix object directly usable by image functions. rawData are stored as well in the map object, for traceability purposes. The neighborhood-related components include the list of neighbor point indices for each kriged data point, as well as the areas of Voronoi polygons associated to all points. Finally the variogram-related components include the VGM model used to simulate the field, as well as the equivalent RandomFields model.

  # Check the mean and standard deviation of generated data
  meanvarSimu(map)

Generate zoning from map for a given probability vector

Given a probability vector, a vector of values is obtained using the quantile function. Display map image with contour levels corresponding to qq values

  qq=quantile(map$krigGrid,na.rm=TRUE,prob=c(0.5,0.7))
  dispZ(map$step,map$krigGrid,valQ=qq)

A zoning is done on the kriged data, by computing the contour lines corresponding to the vector of values given by the probability vector, trimming them to the map boundary and defining zones corresponding to the closed contour lines. A zoning is a list of SpatialPolygons objects. It does not contain data, only polygon geometry.

  # names(ZK):  "resCrit"  "resDist" "resZ" "cL" "qProb"
  ZK=initialZoning(qProb=c(0.5,0.7),map=map)

Plot zoning (14 zones in this case)

   K=ZK$resZ
   Z=K$zonePolygone
   plotZ(Z)

Or a more detailed plot with data underneath - compare it to previous plot of contour lines. We see that contour lines are now extended to the map boundary in order to close zones and that contour lines impossible to close or zones with 0 or 1 point are removed.

   # Outline boundaries of zone 12 with different colors
   dispZ(map$step,map$krigGrid,zonePolygone=Z,iZ=12)
   # Outline all polygons of zone 13 with a different color
   linesSp(Z[[13]],col="blue") # first one in blue
   linesSp(Z[[2]],k=2,col="green") # second one in red, and so on
   # A zone can have one or several holes
   # and each hole is an independent zone. zone 2 has 3 holes (zones 3,8,4). 
   # Due to its shape and to the common borders with the map boundary, 
   # zone 7 is not a hole in zone 2.
   holeSp(Z[[2]])

Junction of 2 zones: Join zone 12 with another zone near by (zone 13). Both zones have the same label

   kmi=optiRG(K,map,12,13,disp=1)
   plotZ(kmi$zonePolygone)

A more detailed plot: each zone is labelled with its number and its mean value

  dispZ(map$step,map$krigGrid,zonePolygone=Z,K=K,boundary=map$boundary,nbLvl=0,id=FALSE,mu=2,cex=0.7) 
  # add quantile values and criterion value for Z
  title(paste(" q=[",toString(round(qq,2)),"]  crit=",round(ZK$resCrit,2),sep="")) 
  # print zone labels
  printLabZ(list(K))
  # print zone areas
  printZsurf(K$zonePolygone)
  # print zone ids
  printZid(Z)
  # remove zones with less than 10 data points
  K=calNei(Z,map$krigData,map$krigSurfVoronoi,map$krigN,nmin=10)
  plotZ(K$zonePolygone)

   val=valZ(map,K)$val
   boxplot(val,col=topo.colors(length(val)))

Generate tree of possible corrections for small zones

2 operations are done for each small zone :

1- remove zone, i.e. merge into englobing zone
2- grow zone

Growing is done in 2 different ways depending on zone proximity to other ones. If zone is isolated (distance to other zones controlled by distIsoZ parameter), it grows bigger but remains isolated from others. Zone growing in that case is performed by finding the contour line close to the current zone contour, that maximizes the zoning quality criterion. A small value of distIsoZ ensures that a small zone have enough space to grow. If zone is non isolated, it is joined to the closest zone with the same label.

  #save all branches resulting from correction steps
  criti<-correctionTree(c(0.4,0.7),map,SAVE=TRUE,ALL=TRUE,LASTPASS=FALSE,distIsoZ=0.01)
  zk=criti$zk

In that case we have for example #4 small zone, hence 3 levels (level 1 is initial zoning, level 2 has 2 branches, level 3 has 4 branches). For each correction step-first branch=zone removal, second-branch=zone junction. The procedure starts with the smallest zone, here zone #4.

   Z21=zk[[2]][[1]]$zonePolygone
   Z22=zk[[2]][[2]]$zonePolygone
   plotZ(Z21,id=TRUE) # result of removal of zone #4
   plotZ(Z22,id=TRUE) # result of growing of zone #4

   # successively:  removal of zone#4  in Z21, growing of zone#4  in Z21
   # removal of zone#4  in Z22, growing of zone#4  in Z22
   # other try with LASTPASS=TRUE removes at last step the zones that are still too small 
   # after all successive corrections
   criti<-correctionTree(c(0.4,0.7),map,SAVE=TRUE,ALL=TRUE,LASTPASS=TRUE,distIsoZ=0.001)

   # other run with ALL=FALSE saves memory by keeping only the first and the last levels
   criti<-correctionTree(c(0.4,0.7),map,SAVE=TRUE,ALL=FALSE,LASTPASS=FALSE,distIsoZ=0.001)