Nothing
R/utilZ.R
In geozoning: Zoning Methods for Spatial Data
##################################################################
#' getZoneId
#'
#' @details get the zone unique identifier
#' @param zone SpatialPolygons
#'
#' @return the zone identifier ( a character vector of length 1)
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' getZoneId(Z[[4]])
#' }
getZoneId=function(zone)
##################################################################
{
id=zone@polygons[[1]]@ID
return(id)
}
###########################################################################
#' Identify
#'
#' @details get the number of a zone with a given identifier in a zoning
#' this is necessary because correction procedures may remove zones from initial#' zoning. Therefore zone numbers change, but identifiers are conserved.
#' @param id zone identifier (character vector)
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return the zone number
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' Identify(6,Z)
#' }
Identify=function(id,Z)
###########################################################################
{
#returns the zone number corresponding to the id given in argument
#returns 0 if no id
for (i in (1:length(Z))){
if (Z[[i]]@polygons[[1]]@ID == id) return(i)
}
return(0)
}
##################################################################
#' getId
#'
#' @details get zone identifier in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#'
#' @return a character vector giving the zone identifier
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' getId(Z,6)
#' }
getId=function(Z,iZ)
##################################################################
{
id = Z[[iZ]]@polygons[[1]]@ID
return(id)
}
##################################################################
#' getIds
#'
#' @details get zone identifiers in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param nums zone numbers
#'
#' @return a character vector giving the zone identifiers
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' geozoning:::getIds(Z)
#' }
getIds=function(Z,nums=NULL)
##################################################################
{
ids=c()
if (is.null(nums))
iZZ = 1:length(Z)
else
iZZ = nums
for (iZ in iZZ)
{
ids = c(ids,Z[[iZ]]@polygons[[1]]@ID)
}
return(ids)
}
##################################################################
#' setId
#'
#' @details assign zone identifier in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#' @param id zone identifier to assign
#'
#' @return a zoning geometry
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' Z1=setId(Z,4,"4")
#' }
setId=function(Z,iZ,id)
##################################################################
{
Z[[iZ]]@polygons[[1]]@ID = as.character(id)
return(Z)
}
##################################################################
#' setIds
#'
#' @details set all zone identifiers in a zoning by assigning zone number to each identifier.
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a zoning geometry
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' Z1=setIds(Z)
#' }
setIds=function(Z)
##################################################################
{
for (iZ in 1:length(Z)){
Z[[iZ]]@polygons[[1]]@ID = as.character(iZ)
}
return(Z)
}
##################################################################
#' maxId
#'
#' @details get highest number corresponding to a zone identifier in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a number
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' geozoning:::maxId(Z)
#' }
maxId=function(Z)
##################################################################
{
m=as.numeric(getId(Z,1))
if (length(Z)<2) return(m)
for (iZ in 2:length(Z)){
id=as.numeric(getId(Z,iZ))
if (id >m) m=id
}
return(as.numeric(m))
}
##################################################################
#' printZid
#'
#' @details print zone identifiers in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return NULL
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' printZid(Z)
printZid = function(Z)
##################################################################
{
for (ii in 1:length(Z))
{
print(paste("ii=",ii," ID=", Z[[ii]]@polygons[[1]]@ID))
}
}
##################################################################
#' getSurf
#'
#' @details get surface area of a given zone within a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#'
#' @return zone area
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' geozoning:::getSurf(Z,1)
#' }
getSurf=function(Z,iZ)
##################################################################
{
return(Z[[iZ]]@polygons[[1]]@Polygons[[1]]@area)
}
##################################################################
#' printZsurf
#'
#' @details print zone surfaces
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param minSize minimum size threshold
#'
#' @return a vector of small zone indices
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' printZsurf(Z,0.03)
printZsurf=function(Z,minSize=0.012)
##################################################################
{
smallZ=c()
gas=c()
for ( iZ in 1:length(Z)){
ga=gArea(Z[[iZ]])
if(ga>=minSize)
print(paste("iZ=",iZ,"area=",round(gArea(Z[[iZ]]),5)))
else{
print(paste("iZ=",iZ," area=",round(gArea(Z[[iZ]]),5)," < minSize(",minSize,")",sep=""))
smallZ=c(smallZ,iZ)
gas=c(gas,ga)
}
}
if (!is.null(gas))
return(smallZ[order(gas)])
else
return(NULL)
}
##################################################################
#' getZonePts
#'
#' @details get all data points within a zone
#' @param ptsp SpatialPointsDataFrame with data values
#' @param zone SpatialPolygons defining a zone
#'
#' @return a list with components
#' \describe{
#' \item{pts}{SpatialPointsDataFrame with the data points within the zone}
#' \item{mask}{Logical vector of the within zone data points indices}
#' }
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' data(mapTest)
#' ptsp=mapTest$krigData
#' res=getZonePts(ptsp,Z[[5]])
#' plotZ(Z)
#' points(res$pts,col="blue",pch=20)
getZonePts=function(ptsp,zone)
##################################################################
{
# pts in zone
IN=numeric()
poly=getPolySp(zone,1)
IN=point.in.polygon(ptsp$x,ptsp$y,poly@coords[,1],poly@coords[,2])
# remove pts in holes
if(nPolySp(zone)>1){
for(k in 2:nPolySp(zone)){
poly=getPolySp(zone,k)
IN=IN-point.in.polygon(ptsp$x,ptsp$y,poly@coords[,1],poly@coords[,2])
}
}
IN=as.logical(IN)
ptsx=ptsp$x[IN]
ptsy=ptsp$y[IN]
ptsd=ptsp[[1]][IN]
ptsub=SpatialPointsDataFrame(coords=cbind(ptsx,ptsy),data=data.frame(z=ptsd))
return(list(pts=ptsub,mask=IN))
}
############################
#' getNumZone
#'
#' @details get zone numbers to which each point in a SpatialPointsDataFrame belongs
#' @param ptsp SpatialPointsDataFrame
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return the zone number
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' getNumZone(mapTest$krigData,Z)
getNumZone=function(ptsp,Z)
############################
{
# ptsp SpatialPointsDataframe
# Z zoning (list of SpatialPolygons)
# ptsp=map$krigData
#
numzone=rep(0,length(ptsp$x))
for (k in 1:length(Z)){
res=getZonePts(ptsp,Z[[k]])
mask=res$mask
numzone[mask]=k
}
return(numzone)
}
##################################################################
#' getNs
#'
#' @details get zone numbers of neighbors of a given zone
#' @param zoneN zone neighborhood Logical matrix
#' @param iZ index of current zone in zoning
#'
#' @return a Logical vector of current zone neighbors
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' K=resZTest
#' Ns=geozoning:::getNs(K$zoneNModif,5) # find neighbors of zone 5
#' }
##################################################################
getNs=function(zoneN,iZ)
{
Ns=zoneN[iZ,]
return(Ns)
}
##################################################################
#' findNptInZone
#'
#' @details find, in a given zone, neighbor points of points belonging to another zone
#' @param K zoning object, as returned by the calNei function
#' @param i1 first zone
#' @param i2 second zone, where to search for neighbors of points in first zone
#' @param map object returned by function genMap
#'
#' @return a two-column matrix, the first column contains indices of pts in first zone which have at least one neighbor in second zone, the second column contains the neighbor indices.
#
#' @export
#'
#' @examples
#' # not run
findNptInZone=function(K,i1,i2,map)
##################################################################
{
ptZ1=K$listZonePoint[[i1]]
ptZ2=K$listZonePoint[[i2]]
ptN=map$krigN
#
# for each pt in Z1
mat=c()
for (ptnum in ptZ1){
# look for neighbor pt in Z2
nei=map$krigN[[ptnum]]
k=match(nei,ptZ2)
n2=ptZ2[k[!is.na(k)]]
if(length(n2)>0){
for (n2k in n2){
mat=rbind(mat,c(ptnum,n2k))
}
}
}
colnames(mat)=c("ptZ1","ptNZ2")
return(mat)
}
##################################################################
#' crComment
#'
#' @details create comment corresponding to holes in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a zoning
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' Z1=geozoning:::crComment(Z)
#' }
crComment = function(Z)
##################################################################
{
# create comments for holes
if (!is.null(Z) & length(Z) >0){
for (iZ in 1:length(Z)){
comment(Z[[iZ]]@polygons[[1]])=createPolygonsComment(Z[[iZ]]@polygons[[1]])
}
}
return(Z)
}
##################################################################
#' testInterSpe
#'
#' @details checks if 2 zones in a zoning share somme common part (using gOverlaps)
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param i1 first zone number
#' @param i2 second zone number
#'
#' @return a Logical value, TRUE if there is an intersection, FALSE if not.
#' @importFrom rgeos gOverlaps
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(Z)
#' Z58=rgeos::gConvexHull(rgeos::gUnion(Z[[8]],Z[[5]]))
#' Z[[length(Z)+1]]=Z58 # add new zone to zoning
#' plotZ(Z)
#' geozoning:::testInterSpe(Z,6,length(Z))
#' }
testInterSpe=function(Z,i1,i2)
##################################################################
# returns TRUE if zones Z[[i1]] and Z[[i2]] share some common part
# but are not within each other
{
#if(gCrosses(Z[[i1]],Z[[i2]])) #problematic sometimes
if(gOverlaps(Z[[i1]],Z[[i2]]))
return(TRUE)
else
return(FALSE)
}
##################################################################
#' testInterSpeZ1
#'
#' @details checks, within a zoning, if a given zone intersects with any other zone not within it
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#'
#' @return a Logical value, TRUE if there is an intersection, FALSE if not.
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(Z)
#' Z58=rgeos::gConvexHull(rgeos::gUnion(Z[[8]],Z[[5]]))
#' Z[[length(Z)+1]]=Z58 # add new zone to zoning
#' plotZ(Z)
#' geozoning:::testInterSpe(Z,6,length(Z))
#' }
testInterSpeZ1=function(Z,iZ)
##################################################################
#returns TRUE if zone Z[[iZ]] intersects with any other zone not within it
{
inter=FALSE
jj=1
while(!inter && (jj<=length(Z)))
{
inter=testInterSpe(Z,iZ,jj)
jj=jj+1
}
return(inter)
}
##################################################################
#' testInterSpeZ
#'
#' @details checks, within a zoning, if any zone intersects with any other zone not within it and not englobing it
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a Logical value, TRUE if there is any intersection, FALSE if not
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(Z)
#' Z58=rgeos::gConvexHull(rgeos::gUnion(Z[[8]],Z[[5]]))
#' Z[[length(Z)+1]]=Z58 # add new zone to zoning
#' plotZ(Z)
#' geozoning:::testInterSpeZ(Z)
#' }
testInterSpeZ=function(Z)
##################################################################
#returns TRUE if any zone intersects with any other zone not within it
{
inter=FALSE
i1=1
while(!inter && (i1<=length(Z)))
{
inter=testInterSpeZ1(Z,i1)
i1=i1+1
}
return(inter)
}
##################################################################
#' printInterZ
#'
#' @details checks intersection of sp and each element of Z
#' @param Z list of zones, each zone is a SpatialPolygons
#' @param sp SpatialPolygons object
#'
#' @return NULL
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' plotZ(Z,mapTest)
#' geozoning:::printInterZ(Z,Z[[1]])
#' }
printInterZ = function(Z,sp)
##################################################################
{
le=length(Z)
if (le == 0) return()
for (i in (1:le))
{
print(paste("i=",i,"intersec=",gIntersects(Z[[i]],sp)))
}
}
##################################################################
#' ptInZone
#'
#' @details description, a paragraph
#' @param zone SpatialPolygons
#' @param pts data points
#' @param numpt data point indices
#'
#' @return 1 if point is within zone, 0 if not
#'
#' @export
#' @importFrom sp point.in.polygon
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' ptInZone(Z[[1]],mapTest$krigData,c(5,500))
ptInZone=function(zone,pts,numpt)
##################################################################
{
spp=SpatialPoints(pts)
logicalPoly=point.in.polygon(spp$x,spp$y,zone@polygons[[1]]@Polygons[[1]]@coords[,1],zone@polygons[[1]]@Polygons[[1]]@coords[,2])
if(length(zone@polygons[[1]]@Polygons)>1) {
#consider pts may be inside holes
for(k in 2:length(zone@polygons[[1]]@Polygons)){
logicalPoly=logicalPoly-point.in.polygon(spp$x,spp$y,zone@polygons[[1]]@Polygons[[k]]@coords[,1],zone@polygons[[1]]@Polygons[[k]]@coords[,2])
}
}
return(logicalPoly[numpt])
}
##################################################################
#' printLabZ
#'
#' @details print zoning labels for a list of zoning objects
#' @param Klist list of zoning objects, typically result of a call to correctionTree
#'
#' @return a list of zoning objects
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' qProb=c(0.1,0.2,0.4);criti=correctionTree(qProb,mapTest) # 2 zonings at last level
#' printLabZ(criti$zk[[2]])
printLabZ=function(Klist)
##################################################################
{
#Klist is a list of several Ks
lk=length(Klist)
labZ=list()
for (k in 1:lk)
{
labk=unlist(Klist[[k]]$lab)
labq=paste(length(unique(labk))-1,"q",sep="")
labZ[[k]]=labk
print(paste(labq,"zone labels=",labZ[k]))
}
return(labZ)
}
##################################################################
#' trLabZone
#'
#' @details transfer zone labels from K1 to K2
#' @param K1 zoning object (such as returned by calNei function)
#' @param K2 zoning object (such as returned by calNei function)
#' @param map object returned by genMap function
#' @param qProb probability vector used to generate quantile values
#' @param disp 0: no info, 1: detailed info
#'
#' @return a zoning object
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' data(resZTest)
#' K=resZTest
#' Ns=geozoning:::getNs(K$zoneNModif,5) # find neighbors of zone 5
#' zf=zoneFusion3(K,5,Ns,mapTest,disp=0) # merge zone 5 with englobing one
#' K2=calNei(zf,mapTest$krigData,mapTest$krigSurfVoronoi,mapTest$krigN)
#' K2=geozoning:::trLabZone(K,K2,mapTest,K$qProb)
#' }
trLabZone=function(K1,K2,map,qProb,disp=0)
##################################################################
# transfer zone labels from K1 to K2
# ids of zones in Z1 and Z2 used to guide transfer
{
Z1=K1$zonePolygone
Z2=K2$zonePolygone
# only K2$lab is modified
# by taking labels from K1
lab2=rep(1,length(Z2))
lab1=K1$lab
id1s=getIds(Z1)
q1 = quantile(map$krigGrid,na.rm=TRUE,prob=qProb)
rate= max(map$krigGrid)-min(map$krigGrid)
for (iZ in 1:length(Z2))
{
id2=getId(Z2,iZ)
numid1=which(id1s==id2)
if (length(numid1) >1)
{
print(id1s)
print(id2)
}
if (length(numid1)!=0)
lab2[iZ]=lab1[numid1]
else #new zone id was created - find its label
{
for (j in 1:length(q1))
{
if (K2$meanZone[iZ]>(q1[j]+0.01*rate))
{
lab2[iZ]=j+1
}
}
}
}
K2$lab=lab2
if(disp)
{
printZid(Z1)
print(K1$lab)
print(K2$lab)
}
return(K2)
}
##################################################################
#' getClosestZone
#'
#' @details get closest non neighbor zone (i.e. excluding neighbor zones and englobing zone)
#' @param iZ current zone number
#' @param Z current zone
#' @param zoneN zone neighborhood Logical matrix
#'
#' @return the closest zone index
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' geozoning:::getClosestZone(4,Z,resZTest$zoneNModif)
#' }
getClosestZone=function(iZ,Z,zoneN)
##################################################################
{
imin=0
ni=1:length(Z)
# exclude current zone
ni=ni[ni != iZ]
# exclude neighbor zones
Ns=getNs(zoneN,iZ)
listeV=grep(TRUE, Ns)
for (i in listeV) ni=ni[ni != i]
# exclude englobing zone
iE = detZoneEng(iZ,Z,zoneN)
#
if(iE != 0) ni = ni[ni !=iE]
# exclude included zones
ir=NULL
for (i in ni)
{
gb = gBuffer(gConvexHull(Z[[iZ]]),byid=TRUE,width=1e-3)
if(gContains(gb,Z[[i]])) ir=c(ir,i)
}
for (i in ir) ni =ni[ni != i]
# compute distances to remaining zones
d0 = 1
for (i in ni)
{
d=gDistance(Z[[iZ]],Z[[i]])
if (d<=d0)
{
imin=i
d0=d
}
}
return(imin)
}
#####################################################################
#' MeanVarWPts
#'
#' @details computes (weighted) mean and variance of zone data
#' @param map object returned by function genMap
#' @param zone SpatialPolygons defining a zone
#' @param w weighting vector (default NULL)
#'
#' @return a list with components
#' \describe{
#' \item{mean}{(weighted) mean of the within zone attribute value}
#' \item{var}{(weighted) variance of the within zone attribute value}
#' }
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' data(mapTest)
#' MeanVarWPts(mapTest,Z[[1]])
#' # Weights are areas of the Voronoi polygons corresponding to data points
#' MeanVarWPts(mapTest,Z[[1]],mapTest$krigSurfVoronoi) #slightly different result
MeanVarWPts=function(map,zone,w=NULL)
#####################################################################
{
ptsp=map$krigData
m=numeric()
v=numeric()
res=getZonePts(ptsp,zone)
mask=res$mask
if (is.null(w))
w=map$krigSurfVoronoi
else
w=rep(1,length(mask))
m=sum(ptsp[[1]][mask]*w[mask]) /sum(w[mask])
d=ptsp[[1]][mask]-m
v=sum(d*d*w[mask]) /sum(w[mask])
return(list(mean=m,var=v))
}
##################################################################
#' nPolyZone
#'
#' @details number of polygons in current zone
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iC current zone number within Z
#'
#' @return the number of polygons within the current zone
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' iC=1;print(paste(nPolyZone(Z,iC),"polygons in zone",iC))
nPolyZone=function(Z,iC)
##################################################################
{
return(length(Z[[iC]]@polygons[[1]]@Polygons))
}
##################################################################
#' maxDistZone
#'
#' @details maximum distance within kth polygon of current zone
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ current zone index
#' @param k polygon number within current zone
#'
#' @return the maximum distance within kth polygon of the current zone
#'
#' @keywords internal
#' @examples
#' \donttest{
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' geozoning:::maxDistZone(Z,5,1)
#' }
maxDistZone=function(Z,iZ,k)
##################################################################
{
return(max(dist(Z[[iZ]]@polygons[[1]]@Polygons[[k]]@coords)))
}
##################################################################
#' getPoly
#'
#' @details get the kth polygon of the current zone in zoning Z
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ current zone index
#' @param k polygon number within current zone
#'
#' @return a polygon (object of class Polygon)
#'
#' @export
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' P1=getPoly(Z,5,1)
#' P2=getPoly(Z,5,2) # second polygon is a hole
#' plot(P1@coords,type="l")
#' lines(P2@coords,type="l",col="blue")
getPoly = function(Z,iZ,k)
##################################################################
{
if (k == 0) return(Z[[iZ]])
# k can be <0 or >0
p=Z[[iZ]]@polygons[[1]]@Polygons[[k]]
return(p)
}
##################################################################
#' polyToSp
#'
#' @details transforms kth polygon of zone into SpatialPolygons
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#' @param k polygon number
#'
#' @return a SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' sph=geozoning:::polyToSp(Z,5,2)
#' plotZ(Z)
#' sp::plot(sph,type="l",col="blue",add=TRUE)
#' }
polyToSp=function(Z,iZ,k)
##################################################################
{
if (k == 0) return(Z[[iZ]])
# k can be <0 or >0
p=Z[[iZ]]@polygons[[1]]@Polygons[[k]]
sp=SpatialPolygons(list(Polygons(list(p),1:1)))
return(sp)
}
##################################################################
#' calcDCrit
#'
#' @details computes distances and criterion value for zoning Z
#' @param Z zoning geometry (list pf SpatialPolygons)
#' @param map object returned by function genMap
#' @param optiCrit criterion choice
#' @param pErr equality tolerance for distance calculations, default 0.9
#' @param simplitol tolerance for spatial polygons geometry simplification, default 0.001
#'
#' @return a list with components
#' \describe{
#' \item{resD}{list with uncorrected and corrected distance matrix}
#' \item{resCrit}{list with criterion and cost values}
#' }
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' Z1=zoneFusion4(Z,6,2)
#' calcDCrit(Z1,mapTest)
calcDCrit=function(Z,map,optiCrit=2,pErr=0.9,simplitol=1e-3)
##################################################################
{
resZ=calNei(Z,map$krigData,map$krigSurfVoronoi,map$krigN,simplitol=simplitol)
le = length(resZ$zonePolygone)
if (le <2) return(list(resD=0,resCrit=0))
resDist=calDistance(typedist=1,map$krigData,resZ$listZonePoint,resZ$zoneN,map$krigSurfVoronoi,resZ$meanZone,pErr=pErr)
resCrit=calCrit(resDist$matDistanceCorr,resZ$zoneNModif,optiCrit)
return(list(resD=resDist,resCrit=resCrit))
}
##################################################################
#' normZcoords
#'
#' @details description, a paragraph
#' @param Z list of SpatialPolygons
#' @param boundary list with components x and y, used to normalize polygons in zoning
#'
#' @return a list with components
#' \describe{
#' \item{Zn}{list of normalized SpatialPolygons}
#' \item{boundaryn}{normalized boundary}
#' }
#' @export
#'
#' @examples
#' # Import shape1 object (was read from a shapefile)
#' shape1 = geozoning::shape1
#' p = shape1@polygons
#' P=sp::SpatialPolygons(p) #SpatialPolygons
#' Z1=list()
#' for (kk in 1:length(P)){Z1[[kk]]=P[kk]} # transform into list of SpatialPolygons
#' bd=list(x=c(7723131,7723132,7723294,7723295,7723131),y=c(3576432,3576814,3576809,3576436,3576432))
#' Z2=normZcoords(Z1,bd)
normZcoords=function(Z,boundary)
##################################################################
{
NZ=length(Z)
Z1=list()
for (iZ in 1:NZ)
{
np=nPolyZone(Z,iZ)
pnl=list()
for (k in 1:np)
{
pk = getPoly(Z,iZ,k)
resn = spnorm(pk,boundary)
if (is.null(resn)) return(NULL)
pnl[[k]] = resn$pn
boundaryn = resn$boundaryn
}
Z1[[iZ]] = SpatialPolygons(list(Polygons(pnl,1:1)))
}
return(list(Zn=Z1,boundaryn=boundaryn))
}
##################################################################
#' createHoles
#'
#' @details description, a paragraph
#' @param Z list of zones, each zone is a SpatialPolygons
#'
#' @return a list of zones where holes are distinct SpatialPolygons
#'
#' @export
#'
#' @examples
#' # not run
createHoles = function(Z)
##################################################################
{
NZ = length(Z)
hole = matrix(NA,nrow=NZ,ncol=NZ)
Z1 = Z
for (iZ in 1:NZ)
{
for (kZ in 1:NZ)
{
if (iZ == kZ) next
h=gContains(Z[[iZ]],Z[[kZ]])
hole[iZ,kZ]=h
if(h)
{
#kz is within iz, create corresp. hole in iz
Z1[[iZ]]=gDifference(Z[[iZ]],Z[[kZ]])
}
} #end for kZ
}
#end for iZ
print(hole)
return(Z1)
}
##################################################################
#' moveHoles
#'
#' @details creates SpatialPolygons excluding holes
#' @param zoneMain SpatialPolygons
#' @param zoneSuppr SpatialPolygons inside main
#'
#' @return a new SpatialPolygons object
#'
#' @keywords internal
#'
#' @examples
#' # not run
moveHoles = function(zoneMain,zoneSuppr)
##################################################################
{
Zone1=zoneMain
le=nPolySp(zoneSuppr)
for ( i in 1:le)
{
poly=zoneSuppr@polygons[[1]]@Polygons[[i]]
polys=Polygons(list(poly),"id")
polySp = SpatialPolygons(list(polys))
if(poly@hole) Zone1=gDifference(Zone1,polySp)
}
return(Zone1)
}
##################################################################
#' find contour for a given quantile value, within an envelope and englobing current zone
#'
#' @details withing a zoning, find contour for a given vRef quantile value, contour contains current zone and is included in envel (spatial Polygon)
#' @param iC zone number
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param K zoning object (such as returned by calNei function)
#' @param map object returned by genMap function
#' @param vRef quantile value
#' @param envel SpatialPolygons within which the contour must be contained
#'
#' @return a list with components
#'\describe{
#' \item{area}{area of SpatialPolygons corresponding to contour}
#' \item{contourSp}{SpatialPolygons corresponding to contour}
#' }
#' @importFrom sp plot
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.3,0.5)
#' criti = correctionTree(qProb,mapTest)
#' best = criti$zk[[2]][[1]]
#' Z=best$zonePolygone
#' plotZ(Z)
#' iC=4
#' envel=calFrame(iC,Z,best$zoneNModif)
#' sp::plot(envel,col="blue",add=TRUE)
#' vRef=quantile(mapTest$krigGrid,0.6)
#' resp=geozoning:::findCinZ(iC,Z,best,mapTest,vRef,envel)
#' sp::plot(resp$contourSp,col="red",add=TRUE)
#' }
findCinZ = function(iC,Z,K,map,vRef,envel)
##################################################################
{
# find contour for a given vRef quantile value
# contour contains Z[[ic]]
# contour included in envel (spatial Polygon)
# init
area=0
#
listContour=list()
listContour = contourAuto(listContour,map$step,map$xsize,map$ysize,map$krigGrid,vRef,map$boundary)
# intersect contour with boundary
# and transform contours into sps
sp=list()
k=0
for (cont in listContour)
{
k=k+1
ps = interCB(cont,map$step,envel=envel)
# returns NULL is contour is degenerate (single point)
if(!is.null(ps)) sp[[k]]=ps else k=k-1
}
if (length(sp)==0) return(NULL)
# check which one contains current zone
spb = sapply(sp,gBuffer,width=1e-3)
gc = sapply(spb,gContains,Z[[iC]])
# test if contour is included in envel
gci = sapply(spb, gWithin, envel)
numc=1:length(sp) # number of non degenerated contours
ga = sapply(sp,gArea)
ind = numc[gc & gci]
# ind may be empty - otherwise take the biggest contour area
if(length(ind)==0) return(NULL)
imax = which(ga[ind]==max(ga[ind]))
imax = imax[1]
im = numc[ind[imax]]
area = max(ga[ind])
contourSp = sp[[im]]
return(list(area=area,contourSp=contourSp))
}
######################
#' modlm
#'
#' @details description, a paragraph
#' @param ptsp SpatialPointsDataFrame with data values
#' @param Z zoning (list of SpatialPolygons)
#'
#' @return the result of a call to lm (anova model with zone number as factor)
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' data(mapTest)
#' ptsp=mapTest$krigData
#' geozoning:::modlm(ptsp,Z)
#' }
modlm=function(ptsp,Z)
######################
{
# ptssp SpatialPointsDataframe
# Z zoning (list of SpatialPolygons)
# ptsp=map$krigData
numz=getNumZone(ptsp,Z)
# add numz to SpatialPointsDataframe
ptsp[[2]]=numz
colnames(ptsp@data)=c("z","numz")
reslm=lm(z~numz,data=ptsp@data)
return(reslm)
}
##################################################################
#' getClosePt
#'
#' @details description, a paragraph
#' @param Z zoning (list of SpatialPolygons)
#' @param iC current zone indes
#' @param iZC close zone index
#' @param disp information level (FALSE-no info)
#'
#' @return a SpatialPoints of length 1
#'
#' @keywords internal
#' @importFrom sp SpatialPoints SpatialPointsDataFrame
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' geozoning:::getClosePt(Z,1,3)
#' plotZ(Z)
#' points( geozoning:::getClosePt(Z,1,3),col="blue",pch=20)
#' }
getClosePt=function(Z,iC,iZC,disp=FALSE)
##################################################################
{
#pt coordinates
a=SpatialPoints(Z[[iC]]@polygons[[1]]@Polygons[[1]]@coords)
b=SpatialPoints(Z[[iZC]]@polygons[[1]]@Polygons[[1]]@coords)
#all distances between pairs of points - lengthy
Fdist= list()
for (i in 1:length(a))
{
Fdist[[i]]<-gDistance(a[i,],b,byid=TRUE)
}
#get the point in zone iZC which is the closest to zone iC
min.dist <- unlist(lapply(Fdist, FUN=function(x) which(x == min(x))[1]))
PolyDist <- unlist(lapply(Fdist, FUN=function(x) min(x)[1]))
pZC = min.dist[which.min(PolyDist)]
if (disp)
{
print(b[pZC])
}
return (b[pZC])
}
########################
#' findZCenter
#'
#' @details find point within zone for pretty labelling
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param num zone number
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a SpatialPoints
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7,
#' # saving initial zoning and last level zonings
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygone
#' findZCenter(Z)
findZCenter=function(Z,num=NULL)
################################
{
if(is.null(num)) num=1:length(Z)
ptz=NULL
for (jj in num){
erosion = Z[[jj]]
width = 0
while(gArea(erosion)>10^-3){
width = width + 0.001
erosion = gBuffer(Z[[jj]],width =-width)
}
ptz=rbind(ptz,coordinates(erosion))
}
rownames(ptz)=NULL
colnames(ptz)=c("x","y")
return(ptz)
}
########################
#' findZCenterpt
#'
#' @details find point within zone for pretty labelling
#' @param data SpatialPointsDataFrame
#' @param K zoning object, as returned by the calNei function
#' @param num zone number or NULL for all zones
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a matrix of x and y coordinates for chosen points with as many rows as zones
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7,
#' # saving initial zoning and last level zonings
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=TRUE)
#' K=criti$zk[[2]][[1]]
#' data=mapTest$krigData
#' geozoning:::findZCenterpt(data,K)
#' }
findZCenterpt=function(data,K,num=NULL)
#######################################
{
Z=K$zonePolygone
if(is.null(num)) num=1:length(Z)
ptz=NULL
for (jj in num){
ipt=K$listZonePoint[[jj]] # indices of within zone pts
ptc=data[ipt,]
gd=rep(NA,length(ptc))
for (k in 1:length(ptc)) {gd[k]=gDistance(ptc[k,],gBoundary(Z[[jj]]))} #smallest distance from each point to zone boundary
ind=which(gd==max(gd)) # maximum distance to zone boundary
ind=ind[1] # if ties
ptz=rbind(ptz,coordinates(ptc[ind,]))
}
rownames(ptz)=NULL
colnames(ptz)=c("x","y")
return(ptz)
}
########################
#' getZsize
#'
#' @details compute maximum x and y values of zoning Z
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a vector with x and y maximum values
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' geozoning:::getZsize(Z)
#' }
getZsize=function(Z)
########################
{
ms=sapply(Z,function(x)return(x@bbox[,"max"]))
vs=apply(ms,1,max)
return(vs)
}
##################################################################
#' interZoneC
#'
#' @details finds two intersection points of a circle with a zone. The circle radius is chosen so that it will intersect both zones given as arguments.
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iC zone number
#' @param iZC other zone number
#' @param closePt SpatialPoints object in other zone used as circle center
#'
#' @return a list with components
#' \describe{
#' \item{spi}{Two SpatialPoints to be used for the junction of the two zones}
#' \item{ord}{Order in which to use the points}
#' }
#' @importFrom rgeos gBoundary gCrosses gEnvelope gIntersection gLength
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(K$zonePolygone) # zoning
#' closePt = geozoning:::getClosePt(Z,6,8)
#' points(closePt,col="red")
#' res = geozoning:::interZoneC(Z,6,8,closePt)
#' points(res$spi,col="red")
#' }
interZoneC=function(Z,iC, iZC, closePt)
##################################################################
{
# returns spatialpoints = 2 points for junction of zones iC and iZC
# iC (current zone) and iZC (close zone)
if (is.null(closePt)) return(NULL) #protection
# coordinates of close zone
b2=gBoundary(Z[[iZC]])
co1 = Z[[iC]]@polygons[[1]]@Polygons[[1]]@coords
co2sp = SpatialPoints(Z[[iZC]]@polygons[[1]]@Polygons[[1]]@coords)
end = TRUE
# draw circle with closePt as center and radius chosen so that it will intersect
x = closePt$x
y = closePt$y
n = 100
r= 0.75*max (dist(co1))
pts = seq(0, 2 * pi, length.out = n)
xy = cbind(x + r * sin(pts), y + r * cos(pts))
circle= SpatialLines(list(Lines(list(Line(xy)), "line")))
#lines(circle)
# intersect circle with close zone
pInter = gIntersection(circle,b2)
spi = NULL
ord=NULL
if (!is.null(pInter) & length(pInter)>=2){
# find intersection pts
i1=which.min(gDistance(pInter[1],co2sp,byid=TRUE))
i2=which.min(gDistance(pInter[2],co2sp,byid=TRUE))
#transform into Spatial Points
spi= co2sp[c(i1,i2),]
# order of points
if (i1 <i2){
ord = (i1+1):(i2-1)
end = FALSE
} else {
ord = (i2+1):(i1-1)
}
nbcoord = length(co2sp)
if (length(ord) >nbcoord/2){
if (!end) {
ord = c((i2+1):(nbcoord),1:(i1-1))
} else {
ord = c((i1+1):(nbcoord),1:(i2-1))
}
}
}
return(list(spi=spi,ord=ord))
}
Try the geozoning package in your browser
Any scripts or data that you put into this service are public.
geozoning documentation built on May 2, 2019, 9:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##################################################################
#' getZoneId
#'
#' @details get the zone unique identifier
#' @param zone SpatialPolygons
#'
#' @return the zone identifier ( a character vector of length 1)
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' getZoneId(Z[[4]])
#' }
getZoneId=function(zone)
##################################################################
{
id=zone@polygons[[1]]@ID
return(id)
}
###########################################################################
#' Identify
#'
#' @details get the number of a zone with a given identifier in a zoning
#' this is necessary because correction procedures may remove zones from initial#' zoning. Therefore zone numbers change, but identifiers are conserved.
#' @param id zone identifier (character vector)
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return the zone number
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' Identify(6,Z)
#' }
Identify=function(id,Z)
###########################################################################
{
#returns the zone number corresponding to the id given in argument
#returns 0 if no id
for (i in (1:length(Z))){
if (Z[[i]]@polygons[[1]]@ID == id) return(i)
}
return(0)
}
##################################################################
#' getId
#'
#' @details get zone identifier in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#'
#' @return a character vector giving the zone identifier
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' getId(Z,6)
#' }
getId=function(Z,iZ)
##################################################################
{
id = Z[[iZ]]@polygons[[1]]@ID
return(id)
}
##################################################################
#' getIds
#'
#' @details get zone identifiers in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param nums zone numbers
#'
#' @return a character vector giving the zone identifiers
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' geozoning:::getIds(Z)
#' }
getIds=function(Z,nums=NULL)
##################################################################
{
ids=c()
if (is.null(nums))
iZZ = 1:length(Z)
else
iZZ = nums
for (iZ in iZZ)
{
ids = c(ids,Z[[iZ]]@polygons[[1]]@ID)
}
return(ids)
}
##################################################################
#' setId
#'
#' @details assign zone identifier in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#' @param id zone identifier to assign
#'
#' @return a zoning geometry
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' Z1=setId(Z,4,"4")
#' }
setId=function(Z,iZ,id)
##################################################################
{
Z[[iZ]]@polygons[[1]]@ID = as.character(id)
return(Z)
}
##################################################################
#' setIds
#'
#' @details set all zone identifiers in a zoning by assigning zone number to each identifier.
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a zoning geometry
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' Z1=setIds(Z)
#' }
setIds=function(Z)
##################################################################
{
for (iZ in 1:length(Z)){
Z[[iZ]]@polygons[[1]]@ID = as.character(iZ)
}
return(Z)
}
##################################################################
#' maxId
#'
#' @details get highest number corresponding to a zone identifier in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a number
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' criti=correctionTree(c(0.4,0.5),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygon
#' geozoning:::maxId(Z)
#' }
maxId=function(Z)
##################################################################
{
m=as.numeric(getId(Z,1))
if (length(Z)<2) return(m)
for (iZ in 2:length(Z)){
id=as.numeric(getId(Z,iZ))
if (id >m) m=id
}
return(as.numeric(m))
}
##################################################################
#' printZid
#'
#' @details print zone identifiers in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return NULL
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' printZid(Z)
printZid = function(Z)
##################################################################
{
for (ii in 1:length(Z))
{
print(paste("ii=",ii," ID=", Z[[ii]]@polygons[[1]]@ID))
}
}
##################################################################
#' getSurf
#'
#' @details get surface area of a given zone within a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#'
#' @return zone area
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' geozoning:::getSurf(Z,1)
#' }
getSurf=function(Z,iZ)
##################################################################
{
return(Z[[iZ]]@polygons[[1]]@Polygons[[1]]@area)
}
##################################################################
#' printZsurf
#'
#' @details print zone surfaces
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param minSize minimum size threshold
#'
#' @return a vector of small zone indices
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' printZsurf(Z,0.03)
printZsurf=function(Z,minSize=0.012)
##################################################################
{
smallZ=c()
gas=c()
for ( iZ in 1:length(Z)){
ga=gArea(Z[[iZ]])
if(ga>=minSize)
print(paste("iZ=",iZ,"area=",round(gArea(Z[[iZ]]),5)))
else{
print(paste("iZ=",iZ," area=",round(gArea(Z[[iZ]]),5)," < minSize(",minSize,")",sep=""))
smallZ=c(smallZ,iZ)
gas=c(gas,ga)
}
}
if (!is.null(gas))
return(smallZ[order(gas)])
else
return(NULL)
}
##################################################################
#' getZonePts
#'
#' @details get all data points within a zone
#' @param ptsp SpatialPointsDataFrame with data values
#' @param zone SpatialPolygons defining a zone
#'
#' @return a list with components
#' \describe{
#' \item{pts}{SpatialPointsDataFrame with the data points within the zone}
#' \item{mask}{Logical vector of the within zone data points indices}
#' }
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' data(mapTest)
#' ptsp=mapTest$krigData
#' res=getZonePts(ptsp,Z[[5]])
#' plotZ(Z)
#' points(res$pts,col="blue",pch=20)
getZonePts=function(ptsp,zone)
##################################################################
{
# pts in zone
IN=numeric()
poly=getPolySp(zone,1)
IN=point.in.polygon(ptsp$x,ptsp$y,poly@coords[,1],poly@coords[,2])
# remove pts in holes
if(nPolySp(zone)>1){
for(k in 2:nPolySp(zone)){
poly=getPolySp(zone,k)
IN=IN-point.in.polygon(ptsp$x,ptsp$y,poly@coords[,1],poly@coords[,2])
}
}
IN=as.logical(IN)
ptsx=ptsp$x[IN]
ptsy=ptsp$y[IN]
ptsd=ptsp[[1]][IN]
ptsub=SpatialPointsDataFrame(coords=cbind(ptsx,ptsy),data=data.frame(z=ptsd))
return(list(pts=ptsub,mask=IN))
}
############################
#' getNumZone
#'
#' @details get zone numbers to which each point in a SpatialPointsDataFrame belongs
#' @param ptsp SpatialPointsDataFrame
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return the zone number
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' getNumZone(mapTest$krigData,Z)
getNumZone=function(ptsp,Z)
############################
{
# ptsp SpatialPointsDataframe
# Z zoning (list of SpatialPolygons)
# ptsp=map$krigData
#
numzone=rep(0,length(ptsp$x))
for (k in 1:length(Z)){
res=getZonePts(ptsp,Z[[k]])
mask=res$mask
numzone[mask]=k
}
return(numzone)
}
##################################################################
#' getNs
#'
#' @details get zone numbers of neighbors of a given zone
#' @param zoneN zone neighborhood Logical matrix
#' @param iZ index of current zone in zoning
#'
#' @return a Logical vector of current zone neighbors
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' K=resZTest
#' Ns=geozoning:::getNs(K$zoneNModif,5) # find neighbors of zone 5
#' }
##################################################################
getNs=function(zoneN,iZ)
{
Ns=zoneN[iZ,]
return(Ns)
}
##################################################################
#' findNptInZone
#'
#' @details find, in a given zone, neighbor points of points belonging to another zone
#' @param K zoning object, as returned by the calNei function
#' @param i1 first zone
#' @param i2 second zone, where to search for neighbors of points in first zone
#' @param map object returned by function genMap
#'
#' @return a two-column matrix, the first column contains indices of pts in first zone which have at least one neighbor in second zone, the second column contains the neighbor indices.
#
#' @export
#'
#' @examples
#' # not run
findNptInZone=function(K,i1,i2,map)
##################################################################
{
ptZ1=K$listZonePoint[[i1]]
ptZ2=K$listZonePoint[[i2]]
ptN=map$krigN
#
# for each pt in Z1
mat=c()
for (ptnum in ptZ1){
# look for neighbor pt in Z2
nei=map$krigN[[ptnum]]
k=match(nei,ptZ2)
n2=ptZ2[k[!is.na(k)]]
if(length(n2)>0){
for (n2k in n2){
mat=rbind(mat,c(ptnum,n2k))
}
}
}
colnames(mat)=c("ptZ1","ptNZ2")
return(mat)
}
##################################################################
#' crComment
#'
#' @details create comment corresponding to holes in a zoning
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a zoning
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' Z1=geozoning:::crComment(Z)
#' }
crComment = function(Z)
##################################################################
{
# create comments for holes
if (!is.null(Z) & length(Z) >0){
for (iZ in 1:length(Z)){
comment(Z[[iZ]]@polygons[[1]])=createPolygonsComment(Z[[iZ]]@polygons[[1]])
}
}
return(Z)
}
##################################################################
#' testInterSpe
#'
#' @details checks if 2 zones in a zoning share somme common part (using gOverlaps)
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param i1 first zone number
#' @param i2 second zone number
#'
#' @return a Logical value, TRUE if there is an intersection, FALSE if not.
#' @importFrom rgeos gOverlaps
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(Z)
#' Z58=rgeos::gConvexHull(rgeos::gUnion(Z[[8]],Z[[5]]))
#' Z[[length(Z)+1]]=Z58 # add new zone to zoning
#' plotZ(Z)
#' geozoning:::testInterSpe(Z,6,length(Z))
#' }
testInterSpe=function(Z,i1,i2)
##################################################################
# returns TRUE if zones Z[[i1]] and Z[[i2]] share some common part
# but are not within each other
{
#if(gCrosses(Z[[i1]],Z[[i2]])) #problematic sometimes
if(gOverlaps(Z[[i1]],Z[[i2]]))
return(TRUE)
else
return(FALSE)
}
##################################################################
#' testInterSpeZ1
#'
#' @details checks, within a zoning, if a given zone intersects with any other zone not within it
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#'
#' @return a Logical value, TRUE if there is an intersection, FALSE if not.
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(Z)
#' Z58=rgeos::gConvexHull(rgeos::gUnion(Z[[8]],Z[[5]]))
#' Z[[length(Z)+1]]=Z58 # add new zone to zoning
#' plotZ(Z)
#' geozoning:::testInterSpe(Z,6,length(Z))
#' }
testInterSpeZ1=function(Z,iZ)
##################################################################
#returns TRUE if zone Z[[iZ]] intersects with any other zone not within it
{
inter=FALSE
jj=1
while(!inter && (jj<=length(Z)))
{
inter=testInterSpe(Z,iZ,jj)
jj=jj+1
}
return(inter)
}
##################################################################
#' testInterSpeZ
#'
#' @details checks, within a zoning, if any zone intersects with any other zone not within it and not englobing it
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @return a Logical value, TRUE if there is any intersection, FALSE if not
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(Z)
#' Z58=rgeos::gConvexHull(rgeos::gUnion(Z[[8]],Z[[5]]))
#' Z[[length(Z)+1]]=Z58 # add new zone to zoning
#' plotZ(Z)
#' geozoning:::testInterSpeZ(Z)
#' }
testInterSpeZ=function(Z)
##################################################################
#returns TRUE if any zone intersects with any other zone not within it
{
inter=FALSE
i1=1
while(!inter && (i1<=length(Z)))
{
inter=testInterSpeZ1(Z,i1)
i1=i1+1
}
return(inter)
}
##################################################################
#' printInterZ
#'
#' @details checks intersection of sp and each element of Z
#' @param Z list of zones, each zone is a SpatialPolygons
#' @param sp SpatialPolygons object
#'
#' @return NULL
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' plotZ(Z,mapTest)
#' geozoning:::printInterZ(Z,Z[[1]])
#' }
printInterZ = function(Z,sp)
##################################################################
{
le=length(Z)
if (le == 0) return()
for (i in (1:le))
{
print(paste("i=",i,"intersec=",gIntersects(Z[[i]],sp)))
}
}
##################################################################
#' ptInZone
#'
#' @details description, a paragraph
#' @param zone SpatialPolygons
#' @param pts data points
#' @param numpt data point indices
#'
#' @return 1 if point is within zone, 0 if not
#'
#' @export
#' @importFrom sp point.in.polygon
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' ptInZone(Z[[1]],mapTest$krigData,c(5,500))
ptInZone=function(zone,pts,numpt)
##################################################################
{
spp=SpatialPoints(pts)
logicalPoly=point.in.polygon(spp$x,spp$y,zone@polygons[[1]]@Polygons[[1]]@coords[,1],zone@polygons[[1]]@Polygons[[1]]@coords[,2])
if(length(zone@polygons[[1]]@Polygons)>1) {
#consider pts may be inside holes
for(k in 2:length(zone@polygons[[1]]@Polygons)){
logicalPoly=logicalPoly-point.in.polygon(spp$x,spp$y,zone@polygons[[1]]@Polygons[[k]]@coords[,1],zone@polygons[[1]]@Polygons[[k]]@coords[,2])
}
}
return(logicalPoly[numpt])
}
##################################################################
#' printLabZ
#'
#' @details print zoning labels for a list of zoning objects
#' @param Klist list of zoning objects, typically result of a call to correctionTree
#'
#' @return a list of zoning objects
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' qProb=c(0.1,0.2,0.4);criti=correctionTree(qProb,mapTest) # 2 zonings at last level
#' printLabZ(criti$zk[[2]])
printLabZ=function(Klist)
##################################################################
{
#Klist is a list of several Ks
lk=length(Klist)
labZ=list()
for (k in 1:lk)
{
labk=unlist(Klist[[k]]$lab)
labq=paste(length(unique(labk))-1,"q",sep="")
labZ[[k]]=labk
print(paste(labq,"zone labels=",labZ[k]))
}
return(labZ)
}
##################################################################
#' trLabZone
#'
#' @details transfer zone labels from K1 to K2
#' @param K1 zoning object (such as returned by calNei function)
#' @param K2 zoning object (such as returned by calNei function)
#' @param map object returned by genMap function
#' @param qProb probability vector used to generate quantile values
#' @param disp 0: no info, 1: detailed info
#'
#' @return a zoning object
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' data(resZTest)
#' K=resZTest
#' Ns=geozoning:::getNs(K$zoneNModif,5) # find neighbors of zone 5
#' zf=zoneFusion3(K,5,Ns,mapTest,disp=0) # merge zone 5 with englobing one
#' K2=calNei(zf,mapTest$krigData,mapTest$krigSurfVoronoi,mapTest$krigN)
#' K2=geozoning:::trLabZone(K,K2,mapTest,K$qProb)
#' }
trLabZone=function(K1,K2,map,qProb,disp=0)
##################################################################
# transfer zone labels from K1 to K2
# ids of zones in Z1 and Z2 used to guide transfer
{
Z1=K1$zonePolygone
Z2=K2$zonePolygone
# only K2$lab is modified
# by taking labels from K1
lab2=rep(1,length(Z2))
lab1=K1$lab
id1s=getIds(Z1)
q1 = quantile(map$krigGrid,na.rm=TRUE,prob=qProb)
rate= max(map$krigGrid)-min(map$krigGrid)
for (iZ in 1:length(Z2))
{
id2=getId(Z2,iZ)
numid1=which(id1s==id2)
if (length(numid1) >1)
{
print(id1s)
print(id2)
}
if (length(numid1)!=0)
lab2[iZ]=lab1[numid1]
else #new zone id was created - find its label
{
for (j in 1:length(q1))
{
if (K2$meanZone[iZ]>(q1[j]+0.01*rate))
{
lab2[iZ]=j+1
}
}
}
}
K2$lab=lab2
if(disp)
{
printZid(Z1)
print(K1$lab)
print(K2$lab)
}
return(K2)
}
##################################################################
#' getClosestZone
#'
#' @details get closest non neighbor zone (i.e. excluding neighbor zones and englobing zone)
#' @param iZ current zone number
#' @param Z current zone
#' @param zoneN zone neighborhood Logical matrix
#'
#' @return the closest zone index
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' geozoning:::getClosestZone(4,Z,resZTest$zoneNModif)
#' }
getClosestZone=function(iZ,Z,zoneN)
##################################################################
{
imin=0
ni=1:length(Z)
# exclude current zone
ni=ni[ni != iZ]
# exclude neighbor zones
Ns=getNs(zoneN,iZ)
listeV=grep(TRUE, Ns)
for (i in listeV) ni=ni[ni != i]
# exclude englobing zone
iE = detZoneEng(iZ,Z,zoneN)
#
if(iE != 0) ni = ni[ni !=iE]
# exclude included zones
ir=NULL
for (i in ni)
{
gb = gBuffer(gConvexHull(Z[[iZ]]),byid=TRUE,width=1e-3)
if(gContains(gb,Z[[i]])) ir=c(ir,i)
}
for (i in ir) ni =ni[ni != i]
# compute distances to remaining zones
d0 = 1
for (i in ni)
{
d=gDistance(Z[[iZ]],Z[[i]])
if (d<=d0)
{
imin=i
d0=d
}
}
return(imin)
}
#####################################################################
#' MeanVarWPts
#'
#' @details computes (weighted) mean and variance of zone data
#' @param map object returned by function genMap
#' @param zone SpatialPolygons defining a zone
#' @param w weighting vector (default NULL)
#'
#' @return a list with components
#' \describe{
#' \item{mean}{(weighted) mean of the within zone attribute value}
#' \item{var}{(weighted) variance of the within zone attribute value}
#' }
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' data(mapTest)
#' MeanVarWPts(mapTest,Z[[1]])
#' # Weights are areas of the Voronoi polygons corresponding to data points
#' MeanVarWPts(mapTest,Z[[1]],mapTest$krigSurfVoronoi) #slightly different result
MeanVarWPts=function(map,zone,w=NULL)
#####################################################################
{
ptsp=map$krigData
m=numeric()
v=numeric()
res=getZonePts(ptsp,zone)
mask=res$mask
if (is.null(w))
w=map$krigSurfVoronoi
else
w=rep(1,length(mask))
m=sum(ptsp[[1]][mask]*w[mask]) /sum(w[mask])
d=ptsp[[1]][mask]-m
v=sum(d*d*w[mask]) /sum(w[mask])
return(list(mean=m,var=v))
}
##################################################################
#' nPolyZone
#'
#' @details number of polygons in current zone
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iC current zone number within Z
#'
#' @return the number of polygons within the current zone
#'
#' @export
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' iC=1;print(paste(nPolyZone(Z,iC),"polygons in zone",iC))
nPolyZone=function(Z,iC)
##################################################################
{
return(length(Z[[iC]]@polygons[[1]]@Polygons))
}
##################################################################
#' maxDistZone
#'
#' @details maximum distance within kth polygon of current zone
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ current zone index
#' @param k polygon number within current zone
#'
#' @return the maximum distance within kth polygon of the current zone
#'
#' @keywords internal
#' @examples
#' \donttest{
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' geozoning:::maxDistZone(Z,5,1)
#' }
maxDistZone=function(Z,iZ,k)
##################################################################
{
return(max(dist(Z[[iZ]]@polygons[[1]]@Polygons[[k]]@coords)))
}
##################################################################
#' getPoly
#'
#' @details get the kth polygon of the current zone in zoning Z
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ current zone index
#' @param k polygon number within current zone
#'
#' @return a polygon (object of class Polygon)
#'
#' @export
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' P1=getPoly(Z,5,1)
#' P2=getPoly(Z,5,2) # second polygon is a hole
#' plot(P1@coords,type="l")
#' lines(P2@coords,type="l",col="blue")
getPoly = function(Z,iZ,k)
##################################################################
{
if (k == 0) return(Z[[iZ]])
# k can be <0 or >0
p=Z[[iZ]]@polygons[[1]]@Polygons[[k]]
return(p)
}
##################################################################
#' polyToSp
#'
#' @details transforms kth polygon of zone into SpatialPolygons
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iZ zone number
#' @param k polygon number
#'
#' @return a SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' sph=geozoning:::polyToSp(Z,5,2)
#' plotZ(Z)
#' sp::plot(sph,type="l",col="blue",add=TRUE)
#' }
polyToSp=function(Z,iZ,k)
##################################################################
{
if (k == 0) return(Z[[iZ]])
# k can be <0 or >0
p=Z[[iZ]]@polygons[[1]]@Polygons[[k]]
sp=SpatialPolygons(list(Polygons(list(p),1:1)))
return(sp)
}
##################################################################
#' calcDCrit
#'
#' @details computes distances and criterion value for zoning Z
#' @param Z zoning geometry (list pf SpatialPolygons)
#' @param map object returned by function genMap
#' @param optiCrit criterion choice
#' @param pErr equality tolerance for distance calculations, default 0.9
#' @param simplitol tolerance for spatial polygons geometry simplification, default 0.001
#'
#' @return a list with components
#' \describe{
#' \item{resD}{list with uncorrected and corrected distance matrix}
#' \item{resCrit}{list with criterion and cost values}
#' }
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' Z1=zoneFusion4(Z,6,2)
#' calcDCrit(Z1,mapTest)
calcDCrit=function(Z,map,optiCrit=2,pErr=0.9,simplitol=1e-3)
##################################################################
{
resZ=calNei(Z,map$krigData,map$krigSurfVoronoi,map$krigN,simplitol=simplitol)
le = length(resZ$zonePolygone)
if (le <2) return(list(resD=0,resCrit=0))
resDist=calDistance(typedist=1,map$krigData,resZ$listZonePoint,resZ$zoneN,map$krigSurfVoronoi,resZ$meanZone,pErr=pErr)
resCrit=calCrit(resDist$matDistanceCorr,resZ$zoneNModif,optiCrit)
return(list(resD=resDist,resCrit=resCrit))
}
##################################################################
#' normZcoords
#'
#' @details description, a paragraph
#' @param Z list of SpatialPolygons
#' @param boundary list with components x and y, used to normalize polygons in zoning
#'
#' @return a list with components
#' \describe{
#' \item{Zn}{list of normalized SpatialPolygons}
#' \item{boundaryn}{normalized boundary}
#' }
#' @export
#'
#' @examples
#' # Import shape1 object (was read from a shapefile)
#' shape1 = geozoning::shape1
#' p = shape1@polygons
#' P=sp::SpatialPolygons(p) #SpatialPolygons
#' Z1=list()
#' for (kk in 1:length(P)){Z1[[kk]]=P[kk]} # transform into list of SpatialPolygons
#' bd=list(x=c(7723131,7723132,7723294,7723295,7723131),y=c(3576432,3576814,3576809,3576436,3576432))
#' Z2=normZcoords(Z1,bd)
normZcoords=function(Z,boundary)
##################################################################
{
NZ=length(Z)
Z1=list()
for (iZ in 1:NZ)
{
np=nPolyZone(Z,iZ)
pnl=list()
for (k in 1:np)
{
pk = getPoly(Z,iZ,k)
resn = spnorm(pk,boundary)
if (is.null(resn)) return(NULL)
pnl[[k]] = resn$pn
boundaryn = resn$boundaryn
}
Z1[[iZ]] = SpatialPolygons(list(Polygons(pnl,1:1)))
}
return(list(Zn=Z1,boundaryn=boundaryn))
}
##################################################################
#' createHoles
#'
#' @details description, a paragraph
#' @param Z list of zones, each zone is a SpatialPolygons
#'
#' @return a list of zones where holes are distinct SpatialPolygons
#'
#' @export
#'
#' @examples
#' # not run
createHoles = function(Z)
##################################################################
{
NZ = length(Z)
hole = matrix(NA,nrow=NZ,ncol=NZ)
Z1 = Z
for (iZ in 1:NZ)
{
for (kZ in 1:NZ)
{
if (iZ == kZ) next
h=gContains(Z[[iZ]],Z[[kZ]])
hole[iZ,kZ]=h
if(h)
{
#kz is within iz, create corresp. hole in iz
Z1[[iZ]]=gDifference(Z[[iZ]],Z[[kZ]])
}
} #end for kZ
}
#end for iZ
print(hole)
return(Z1)
}
##################################################################
#' moveHoles
#'
#' @details creates SpatialPolygons excluding holes
#' @param zoneMain SpatialPolygons
#' @param zoneSuppr SpatialPolygons inside main
#'
#' @return a new SpatialPolygons object
#'
#' @keywords internal
#'
#' @examples
#' # not run
moveHoles = function(zoneMain,zoneSuppr)
##################################################################
{
Zone1=zoneMain
le=nPolySp(zoneSuppr)
for ( i in 1:le)
{
poly=zoneSuppr@polygons[[1]]@Polygons[[i]]
polys=Polygons(list(poly),"id")
polySp = SpatialPolygons(list(polys))
if(poly@hole) Zone1=gDifference(Zone1,polySp)
}
return(Zone1)
}
##################################################################
#' find contour for a given quantile value, within an envelope and englobing current zone
#'
#' @details withing a zoning, find contour for a given vRef quantile value, contour contains current zone and is included in envel (spatial Polygon)
#' @param iC zone number
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param K zoning object (such as returned by calNei function)
#' @param map object returned by genMap function
#' @param vRef quantile value
#' @param envel SpatialPolygons within which the contour must be contained
#'
#' @return a list with components
#'\describe{
#' \item{area}{area of SpatialPolygons corresponding to contour}
#' \item{contourSp}{SpatialPolygons corresponding to contour}
#' }
#' @importFrom sp plot
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.3,0.5)
#' criti = correctionTree(qProb,mapTest)
#' best = criti$zk[[2]][[1]]
#' Z=best$zonePolygone
#' plotZ(Z)
#' iC=4
#' envel=calFrame(iC,Z,best$zoneNModif)
#' sp::plot(envel,col="blue",add=TRUE)
#' vRef=quantile(mapTest$krigGrid,0.6)
#' resp=geozoning:::findCinZ(iC,Z,best,mapTest,vRef,envel)
#' sp::plot(resp$contourSp,col="red",add=TRUE)
#' }
findCinZ = function(iC,Z,K,map,vRef,envel)
##################################################################
{
# find contour for a given vRef quantile value
# contour contains Z[[ic]]
# contour included in envel (spatial Polygon)
# init
area=0
#
listContour=list()
listContour = contourAuto(listContour,map$step,map$xsize,map$ysize,map$krigGrid,vRef,map$boundary)
# intersect contour with boundary
# and transform contours into sps
sp=list()
k=0
for (cont in listContour)
{
k=k+1
ps = interCB(cont,map$step,envel=envel)
# returns NULL is contour is degenerate (single point)
if(!is.null(ps)) sp[[k]]=ps else k=k-1
}
if (length(sp)==0) return(NULL)
# check which one contains current zone
spb = sapply(sp,gBuffer,width=1e-3)
gc = sapply(spb,gContains,Z[[iC]])
# test if contour is included in envel
gci = sapply(spb, gWithin, envel)
numc=1:length(sp) # number of non degenerated contours
ga = sapply(sp,gArea)
ind = numc[gc & gci]
# ind may be empty - otherwise take the biggest contour area
if(length(ind)==0) return(NULL)
imax = which(ga[ind]==max(ga[ind]))
imax = imax[1]
im = numc[ind[imax]]
area = max(ga[ind])
contourSp = sp[[im]]
return(list(area=area,contourSp=contourSp))
}
######################
#' modlm
#'
#' @details description, a paragraph
#' @param ptsp SpatialPointsDataFrame with data values
#' @param Z zoning (list of SpatialPolygons)
#'
#' @return the result of a call to lm (anova model with zone number as factor)
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' data(mapTest)
#' ptsp=mapTest$krigData
#' geozoning:::modlm(ptsp,Z)
#' }
modlm=function(ptsp,Z)
######################
{
# ptssp SpatialPointsDataframe
# Z zoning (list of SpatialPolygons)
# ptsp=map$krigData
numz=getNumZone(ptsp,Z)
# add numz to SpatialPointsDataframe
ptsp[[2]]=numz
colnames(ptsp@data)=c("z","numz")
reslm=lm(z~numz,data=ptsp@data)
return(reslm)
}
##################################################################
#' getClosePt
#'
#' @details description, a paragraph
#' @param Z zoning (list of SpatialPolygons)
#' @param iC current zone indes
#' @param iZC close zone index
#' @param disp information level (FALSE-no info)
#'
#' @return a SpatialPoints of length 1
#'
#' @keywords internal
#' @importFrom sp SpatialPoints SpatialPointsDataFrame
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' geozoning:::getClosePt(Z,1,3)
#' plotZ(Z)
#' points( geozoning:::getClosePt(Z,1,3),col="blue",pch=20)
#' }
getClosePt=function(Z,iC,iZC,disp=FALSE)
##################################################################
{
#pt coordinates
a=SpatialPoints(Z[[iC]]@polygons[[1]]@Polygons[[1]]@coords)
b=SpatialPoints(Z[[iZC]]@polygons[[1]]@Polygons[[1]]@coords)
#all distances between pairs of points - lengthy
Fdist= list()
for (i in 1:length(a))
{
Fdist[[i]]<-gDistance(a[i,],b,byid=TRUE)
}
#get the point in zone iZC which is the closest to zone iC
min.dist <- unlist(lapply(Fdist, FUN=function(x) which(x == min(x))[1]))
PolyDist <- unlist(lapply(Fdist, FUN=function(x) min(x)[1]))
pZC = min.dist[which.min(PolyDist)]
if (disp)
{
print(b[pZC])
}
return (b[pZC])
}
########################
#' findZCenter
#'
#' @details find point within zone for pretty labelling
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param num zone number
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a SpatialPoints
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7,
#' # saving initial zoning and last level zonings
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=TRUE)
#' Z=criti$zk[[2]][[1]]$zonePolygone
#' findZCenter(Z)
findZCenter=function(Z,num=NULL)
################################
{
if(is.null(num)) num=1:length(Z)
ptz=NULL
for (jj in num){
erosion = Z[[jj]]
width = 0
while(gArea(erosion)>10^-3){
width = width + 0.001
erosion = gBuffer(Z[[jj]],width =-width)
}
ptz=rbind(ptz,coordinates(erosion))
}
rownames(ptz)=NULL
colnames(ptz)=c("x","y")
return(ptz)
}
########################
#' findZCenterpt
#'
#' @details find point within zone for pretty labelling
#' @param data SpatialPointsDataFrame
#' @param K zoning object, as returned by the calNei function
#' @param num zone number or NULL for all zones
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a matrix of x and y coordinates for chosen points with as many rows as zones
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7,
#' # saving initial zoning and last level zonings
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=TRUE)
#' K=criti$zk[[2]][[1]]
#' data=mapTest$krigData
#' geozoning:::findZCenterpt(data,K)
#' }
findZCenterpt=function(data,K,num=NULL)
#######################################
{
Z=K$zonePolygone
if(is.null(num)) num=1:length(Z)
ptz=NULL
for (jj in num){
ipt=K$listZonePoint[[jj]] # indices of within zone pts
ptc=data[ipt,]
gd=rep(NA,length(ptc))
for (k in 1:length(ptc)) {gd[k]=gDistance(ptc[k,],gBoundary(Z[[jj]]))} #smallest distance from each point to zone boundary
ind=which(gd==max(gd)) # maximum distance to zone boundary
ind=ind[1] # if ties
ptz=rbind(ptz,coordinates(ptc[ind,]))
}
rownames(ptz)=NULL
colnames(ptz)=c("x","y")
return(ptz)
}
########################
#' getZsize
#'
#' @details compute maximum x and y values of zoning Z
#' @param Z zoning geometry (list of SpatialPolygons)
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a vector with x and y maximum values
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' geozoning:::getZsize(Z)
#' }
getZsize=function(Z)
########################
{
ms=sapply(Z,function(x)return(x@bbox[,"max"]))
vs=apply(ms,1,max)
return(vs)
}
##################################################################
#' interZoneC
#'
#' @details finds two intersection points of a circle with a zone. The circle radius is chosen so that it will intersect both zones given as arguments.
#' @param Z zoning geometry (list of SpatialPolygons)
#' @param iC zone number
#' @param iZC other zone number
#' @param closePt SpatialPoints object in other zone used as circle center
#'
#' @return a list with components
#' \describe{
#' \item{spi}{Two SpatialPoints to be used for the junction of the two zones}
#' \item{ord}{Order in which to use the points}
#' }
#' @importFrom rgeos gBoundary gCrosses gEnvelope gIntersection gLength
#'
#' @keywords internal
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' qProb=c(0.2,0.5)
#' ZK = initialZoning(qProb, mapTest)
#' K=ZK$resZ
#' Z=K$zonePolygone
#' plotZ(K$zonePolygone) # zoning
#' closePt = geozoning:::getClosePt(Z,6,8)
#' points(closePt,col="red")
#' res = geozoning:::interZoneC(Z,6,8,closePt)
#' points(res$spi,col="red")
#' }
interZoneC=function(Z,iC, iZC, closePt)
##################################################################
{
# returns spatialpoints = 2 points for junction of zones iC and iZC
# iC (current zone) and iZC (close zone)
if (is.null(closePt)) return(NULL) #protection
# coordinates of close zone
b2=gBoundary(Z[[iZC]])
co1 = Z[[iC]]@polygons[[1]]@Polygons[[1]]@coords
co2sp = SpatialPoints(Z[[iZC]]@polygons[[1]]@Polygons[[1]]@coords)
end = TRUE
# draw circle with closePt as center and radius chosen so that it will intersect
x = closePt$x
y = closePt$y
n = 100
r= 0.75*max (dist(co1))
pts = seq(0, 2 * pi, length.out = n)
xy = cbind(x + r * sin(pts), y + r * cos(pts))
circle= SpatialLines(list(Lines(list(Line(xy)), "line")))
#lines(circle)
# intersect circle with close zone
pInter = gIntersection(circle,b2)
spi = NULL
ord=NULL
if (!is.null(pInter) & length(pInter)>=2){
# find intersection pts
i1=which.min(gDistance(pInter[1],co2sp,byid=TRUE))
i2=which.min(gDistance(pInter[2],co2sp,byid=TRUE))
#transform into Spatial Points
spi= co2sp[c(i1,i2),]
# order of points
if (i1 <i2){
ord = (i1+1):(i2-1)
end = FALSE
} else {
ord = (i2+1):(i1-1)
}
nbcoord = length(co2sp)
if (length(ord) >nbcoord/2){
if (!end) {
ord = c((i2+1):(nbcoord),1:(i1-1))
} else {
ord = c((i1+1):(nbcoord),1:(i2-1))
}
}
}
return(list(spi=spi,ord=ord))
}
Try the geozoning package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.