R/util.R
In hazaeljones/geozoning: Zoning Methods for Spatial Data
###################################################################
#' getCoords
#'
#' @details get SpatialPolygons coordinates
#' @param sp SpatialPolygons
#' @param k polygon number
#'
#' @return some coordinates
#'
#' @keywords internal
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' geozoning:::getCoords(Z[[1]])
getCoords=function(sp,k=1)
##################################################################
{
coords=sp@polygons[[1]]@Polygons[[k]]@coords
return(coords)
}
##################################################################
#' spToSL
#'
#' @details tranform SpatialPolygons into SpatialLines
#' @param sp SpatialPolygons
#'
#' @return a SpatialLines
#'
#' @export
#' @importFrom sp SpatialLines Lines
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' spToSL(Z[[5]])
spToSL=function(sp)
##################################################################
{
co = getCoords(sp)
li = Lines(list(Line(co)),ID="1")
lis = SpatialLines(list(li))
return(lis)
}
##################################################################
#' contourArea
#'
#' @details area corresponding to closed contour line
#' @param co contour line
#'
#' @return the area within the contour line
#'
#' @export
#' @importFrom sp SpatialPolygons SpatialPointsDataFrame Polygons Polygon
#' @importFrom maptools ContourLines2SLDF
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' contourArea(cL[[8]])
contourArea=function(co)
##################################################################
{
contour = ContourLines2SLDF(list(co))
coords = contour@lines[[1]]@Lines[[1]]@coords
poly=Polygon(coords)
polys=Polygons(list(poly),"id")
contourSp = SpatialPolygons(list(polys))
surface = gArea(contourSp)
return(surface)
}
##################################################################
#' listContourArea
#'
#' @details area of all contour lines in list
#' @param cL list of contour lines
#'
#' @return a list of areas
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' geozoning:::listContourArea(cL)
listContourArea=function(cL)
##################################################################
{
return(sapply(cL,contourArea))
}
##################################################################
#' contourToSpp
#'
#' @details transform contour line into SpatialPolygons
#' @param co contour line (list with contour level and x,y coordinates
#' @param step grid resolution
#'
#' @return a list with components
#'\describe{
#'\item{sp}{SpatialPolygons corresponding to contour line}
#'\item{contour}{SpatialLines corresponding to contour line}
#'\item{polyBuff}{SpatialPolygons corresponding to buffer around contour line}
#'\item{surface}{SpatialPolygons area}
#'}
#'
#' @export
#' @importFrom rgeos gBuffer gArea gCentroid gContains gConvexHull gDifference gDistance gIntersects gWithin
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' contourToSpp(cL[[8]],0.1)
contourToSpp=function(co,step)
##################################################################
{
contour = ContourLines2SLDF(list(co))
coords = contour@lines[[1]]@Lines[[1]]@coords
# attention gBuffer renvoie 2 polygones
polyBuff = gBuffer(contour,width=0.0001*(1/step))
poly=Polygon(coords)
polys=Polygons(list(poly),"id")
contourSp = SpatialPolygons(list(polys))
surface = gArea(contourSp)
return(list(sp=contourSp,contour=contour,polyBuff=polyBuff,surface=surface))
}
##################################################################
#' nPolySp
#' @details number of polygons in SpatialPolygons
#' @param sp SpatialPolygons
#' @return the number of polygons within the current zone
#'
#' @keywords internal
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' print(paste(geozoning:::nPolySp(Z[[2]]),"polygons"))
nPolySp =function(sp)
##################################################################
{
return(length(sp@polygons[[1]]@Polygons))
}
##################################################################
#' holeSp
#'
#' @details number of holes in SpatialPolygons
#' @param sp SpatialPolygons
#'
#' @return the number of holes within sp
#'
#' @export
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' holeSp(Z[[5]]) #zone 5 has 1 hole
holeSp = function(sp)
##################################################################
{
le=nPolySp(sp)
nh=0
if (le == 0) return(0)
for (i in 1:le)
{
poly=sp@polygons[[1]]@Polygons[[i]]
if(poly@hole) nh=nh+1
}
return (nh)
}
##################################################################
#' maxDistSP
#'
#' @details maximum distance within kth polygon of current zone
#' @param sp SpatialPolygons
#' @return the maximum distance within sp
#'
#' @keywords internal
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' geozoning:::maxDistSP(Z[[5]])
maxDistSP=function(sp)
##################################################################
{
return(max(dist(sp@polygons[[1]]@Polygons[[1]]@coords)))
}
##################################################################
#' getPolySp
#'
#' @details get the kth polygon of the current SpatialPolygons
#' @param sp SpatialPolygons object
#' @param k polygon number
#'
#' @return a polygon (object of class Polygon)
#'
#' @keywords internal
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' sp=Z[[5]]
#' P1=geozoning:::getPolySp(sp,1)
#' P2=geozoning:::getPolySp(sp,2) # second polygon is a hole
#' plot(P1@coords,type="l")
#' lines(P2@coords,type="l",col="blue")
getPolySp = function(sp,k=1)
##################################################################
{
p=sp@polygons[[1]]@Polygons[[k]]
return(p)
}
##################################################################
#' polyToSp2
#'
#' @details transforms polygon into SpatialPolygons
#' @param p polygon
#'
#' @return a SpatialPolygons
#'
#' @export
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' sp=Z[[5]]
#' P1=geozoning:::getPolySp(sp,1)
#' sph=polyToSp2(P1)
#' plotZ(Z)
#' sp::plot(sph,col="blue",lwd=2,add=TRUE)
polyToSp2=function(p)
##################################################################
{
sp=SpatialPolygons(list(Polygons(list(p),1:1)))
return(sp)
}
##################################################################
#' lineToSp
#'
#' @details transform closed line into SpatialPolygons
#' @param lin list with x and y line coordinates
#'
#' @return a SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' lin=data.frame(x=cL[[8]]$x,y=cL[[8]]$y)
#' sp=geozoning:::lineToSp(lin)
lineToSp=function(lin)
##################################################################
{
sp=SpatialPolygons(list(Polygons(list(Polygon(lin,hole = FALSE)), "1")))
return(sp)
}
##################################################################
#' spnorm
#'
#' @details normalize Polygon according to border limits
#' @param sp object of class Polygons
#' @param boundary list with x and y components, used to normalize sp
#'
#' @return a list with components
#' \describe{
#' \item{pn}{normalized Polygon}
#' \item{boundaryn}{normalized boundary}
#' }
#'
#' @keywords internal
#'
#' @examples
#' # shape1: result of call to readS on shapefile
#' z=geozoning::shape1
#' bb=list(x=z@bbox[1,],y=z@bbox[2,])
#' p=z@polygons
#' p1=p[[1]]
#' P1=p1@Polygons[[1]]
#' NP1=geozoning:::spnorm(P1,bb)$pn
#' Nbb=geozoning:::spnorm(P1,bb)$boundaryn
#' plot(NP1@coords,xlim=Nbb$x,ylim=Nbb$y)
spnorm = function(sp, boundary)
##################################################################
{
xmin=min(boundary$x)
xmax=max(boundary$x)
ymin=min(boundary$y)
ymax=max(boundary$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
if (abs(ymax-ymin) < 1e-6) return(NULL)
x= sp@coords[,1]
y = sp@coords[,2]
x = (x-xmin)/(xmax-xmin)
y = (y-ymin)/(ymax-ymin)
pn = Polygon(cbind(x,y))
bx=boundary$x
by=boundary$y
bx = (bx-xmin)/(xmax-xmin)
by = (by-ymin)/(ymax-ymin)
bn = list(x=bx,y=by)
return(list(pn=pn, boundaryn=bn))
}
##################################################################
#' normalize data coordinates and border
#'
#' @details normalize boundary between 0 and 1 and data coordinates accordingly
#' @param data data frame with x and y components
#' @param bd boundary (list with x and y components)
#'
#' @return a list with components
#'\describe{
#' \item{dataN}{normalized data}
#' \item{boundaryN}{normalized boundary}
#' \item{xmin}{minimum vaue of x within boundary}
#' \item{xmax}{maximum vaue of x within boundary}
#' \item{ymin}{minimum vaue of y within boundary}
#' \item{ymax}{maximum vaue of y within boundary}
#' }
#'
#' @keywords internal
#'
#' @examples
#' x=runif(100, min=0, max=1)
#' y=runif(100, min=0.2, max=1.7)
#' range(x) # not [0,1]
#' tabData=data.frame(x=x,y=y)
#' bd=list(x=c(0,0,1,1,0), y=c(0.2,1.7,1.7,0.2,0.2))
#' res=geozoning:::datanorm(tabData,bd)
#' apply(res$dataN,2,range)#
datanorm = function(data, bd)
##################################################################
{
xmin=min(bd$x)
xmax=max(bd$x)
ymin=min(bd$y)
ymax=max(bd$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
if (abs(ymax-ymin) < 1e-6) return(NULL)
x = data$x
y = data$y
x = (x-xmin)/(xmax-xmin)
y = (y-ymin)/(ymax-ymin)
data$x = x
data$y = y
#normalize border
bx=bd$x
by=bd$y
bx = (bx-xmin)/(xmax-xmin)
by = (by-ymin)/(ymax-ymin)
bdn = list(x=bx,y=by)
return(list(dataN=data, boundaryN=bdn,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax))
}
##################################################################
#' normalize data coords
#'
#' @details normalize data coordinates between 0 and 1 with different ratios for x and y
#' @param data frame with x and y components
#'
#' @return a normalized data frame
#'
#' @keywords internal
#'
#' @examples
#' nPoints=500
#' x=runif(nPoints, min=0, max=1)
#' y=runif(nPoints, min=0, max=1)
#' range(x) # not [0,1]
#' tabData=data.frame(x=x,y=y)
#' tabData=geozoning:::datanormXY(tabData) # x,y ranges are now [0,1]
datanormXY = function(data)
##################################################################
{
xmin=min(data$x)
xmax=max(data$x)
ymin=min(data$y)
ymax=max(data$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
if (abs(ymax-ymin) < 1e-6) return(NULL)
x = data$x
y = data$y
x = (x-xmin)/(xmax-xmin)
y = (y-ymin)/(ymax-ymin)
data$x = x
data$y = y
return(data)
}
##################################################################
#' normalize data coords with same ratio (for non square field)
#'
#' @details normalize x between 0 and 1, y and boundary with same ratio
#' @param data frame with x and y components
#' @param bd list with x and y components
#'
#' @return a list with components
#'\describe{
#' \item{dataN}{normalized data}
#' \item{boundaryN}{normalized boundary}
#' \item{ratio}{normalizing ratio}
#' \item{xmin}{minimum value of x within boundary}
#' \item{xmax}{maximum value of x within boundary}
#' \item{ymin}{minimum value of y within boundary}
#' \item{ymax}{maximum value of y within boundary}
#' }
#'
#' @export
#'
#' @examples
#' x=runif(100, min=0, max=1)
#' y=runif(100, min=0.2, max=1.7)
#' range(x) # not [0,1]
#' tabData=data.frame(x=x,y=y)
#' bd=list(x=c(0,0,1,1,0), y=c(0.2,1.7,1.7,0.2,0.2))
#' res=datanormX(tabData,bd)
#' apply(res$dataN,2,range)# x range is now [0,1], not y range
#' res$ratio # normalization ratio
datanormX = function(data,bd)
##################################################################
{
xmin=min(data$x)
xmax=max(data$x)
ymin=min(data$y)
ymax=max(data$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
x = data$x
y = data$y
ratio=xmax-xmin
x = (x-xmin)/ratio
y = (y-ymin)/ratio
data$x = x
data$y = y
#normalize boarder
bx=bd$x
by=bd$y
bx = (bx-xmin)/ratio
by = (by-ymin)/ratio
bdn = list(x=bx,y=by)
return(list(dataN=data, boundaryN=bdn,ratio=ratio,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax))
}
##################################################################
#' normSize
#'
#' @details normalize thresholds for small zone detection and no grow zone, considering mapo boundary
#' @param boundaryN normalized map boundary
#' @param minSize minimum size threshold
#' @param minSizeNG no grow size threshold
#'
#' @return a list with components
#' \describe{
#' \item{minSize}{normalized minimum size threshold}
#' \item{minSizeNG}{normalized no grow size threshold}
#' }
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' resT=normSize(mapTest$boundary,0.012,0.001)#normalize thresholds relatively to map boundary area
normSize=function(boundaryN,minSize,minSizeNG)
##################################################################
{
# normalize threshold for small zone detection
coords=boundaryN
poly=Polygon(coords)
polys=Polygons(list(poly),"id")
boundarySp = SpatialPolygons(list(polys))
boundaryArea=gArea(boundarySp)
minSize=minSize/boundaryArea
minSizeNG=minSizeNG/boundaryArea
print(paste("after standardization minSize=",minSize,"minSizeNG=",minSizeNG))
return(list(minSize=minSize,minSizeNG=minSizeNG))
}
##################################################################
#' genQseq
#'
#' @details description, a paragraph
#' @param qProb probability vector used to generate quantile values
#' @param K zoning object, as returned by the calNei function
#' @param map object returned by function genMap
#' @param i1 current zone index
#' @param i2 englobing zone index
#' @param LEQ length of quantile sequence
#' @param MAXP maximum shift from center for quantile sequence
#' @param disp 0: no info, 1: some info
#'
#' @return a plot
#'
#' @export
#' @examples
#' data(mapTest)
#' qProb=c(0.4,0.7)
#' ZK=initialZoning(qProb,mapTest)
#' K=ZK$resZ
#' print(K$lab)
#' genQseq(qProb,K,mapTest,1,2) # from label 3 to label 2
genQseq = function(qProb,K,map,i1,i2,LEQ=5,MAXP=0.1,disp=0)
##################################################################
{
# i1 = current zone index, i2= englobing zone index
#
lab1 = K$lab[i1]
lab2 = K$lab[i2]
nq = length(qProb)
inc = TRUE # increase quantile value
if (lab1 >= lab2) inc = FALSE # decrease quantile value
valRef= quantile(map$krigGrid,na.rm=TRUE,prob=qProb)
# find quantile value corresponding to zone i1
q1min = max(1,lab1-1)
q1max = min(lab1,nq)
if(inc)
q1=q1max
else
q1=q1min
prob1 = qProb[q1]
q1= quantile(map$krigGrid,na.rm=TRUE,prob=prob1)
# increase or decrease quantile value
if (inc)
{
q2 = quantile(map$krigGrid,na.rm=TRUE,prob= min(prob1 + MAXP, 0.99))
Qseq =seq(q1,q2,length.out=LEQ+1)
}
else
{
q3 = quantile(map$krigGrid,na.rm=TRUE,prob= max(prob1 - MAXP, 0.01))
Qseq = seq(q1,q3,length.out=LEQ+1)
}
Qseq=Qseq[-1]
Qseq=sort(unique(Qseq))
if (disp>0)
{
print("Qseq=")
print(Qseq)
}
return(Qseq)
}
##################################################################
#' checkContour
#'
#' @details check admissibility for contour line: surface >minSizeNG and refPoint close enough
#' @param contourSp SpatialPolygons corresponding to closed contour line
#' @param step grid resolution
#' @param refPoint referene point
#' @param minSizeNG zone area threshold under which a zone is not admissible
#'
#' @return Null if contour is not admissible or a list with components
#' \describe{
#' \item{contourSp}{SpatialPolygons corresponding to admissible contour}
#' \item{}{polyBuff}{SpatialPolygons corresponding to gBuffer around admissible contour}
#' }
#' @importFrom sp Polygon
#' @importFrom rgeos plot
#' @importFrom rgeos gCentroid
#' @export
#'
#' @examples
#' data(mapTest)
#' cL=contourAuto(list(),mapTest$step,mapTest$xsize,mapTest$ysize,
#' mapTest$krigGrid,c(5,7),mapTest$boundary)
#' pG=polyToSp2(sp::Polygon(mapTest$boundary)) #SpatialPolygons corresponding to map boundary
#' rgeos::plot(pG)
#' sp8 = contourToSpp(cL[[8]],0.1)$sp
#' refPoint = rgeos::gCentroid(sp8)
#' resp=checkContour(sp8,mapTest$step,refPoint)
#' rgeos::plot(resp$contourSp,col="red",add=TRUE)
checkContour = function(contourSp,step,refPoint,minSizeNG=1e-3)
##################################################################
{
# contourSp spatial object
polyBuff = gBuffer(contourSp,width=0.0001*(1/step)) #spatialPolygon
surface = gArea(contourSp)
condi = surface <=minSizeNG || gDistance(gCentroid(polyBuff),refPoint)>=0.1
if (condi)
return (NULL)
else
return (list(contourSp=contourSp,polyBuff=polyBuff))
}
##################################################################
#' cleanSp
#'
#' @details removes from sp polygons that are too small (artefacts of gDifference)
#' @param sp SpatialPolygons
#' @param tol minimum area for removal
#' @return a SpatialPolygons
#'
#' @export
#'
#' @examples
#' # not run
cleanSp = function(sp,tol=1e-5)
##################################################################
{
# number of polygons
n=nPolySp(sp)
polyl=list()
k=0
# remove polygons that are too small (<1e-6)
# (artefacts of gDifference)
for (i in 1:n)
{
poly=sp@polygons[[1]]@Polygons[[i]]
area=poly@area
if (area >= tol)
{
k=k+1
polyl[[k]]=poly
}
}
if (length(polyl)==0) return(NULL)
polys=Polygons(polyl,"id")
spc= SpatialPolygons(list(polys))
return(spc)
}
##################################################################
#' ptInSp
#'
#' @details finds data points in sp
#' @param sp SpatialPolygons
#' @param pts data points
#' @param hole if TRUE also consider points in holes
#'
#' @return a data frame with data points within sp
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' ptsInSp(Z[[5]],mapTest$krigData) # 5 data points within zone 5
ptsInSp=function(sp,pts,hole=FALSE)
##################################################################
{
PointsSpatiaux=SpatialPoints(pts)
logicalPoly=point.in.polygon(PointsSpatiaux$x,PointsSpatiaux$y,sp@polygons[[1]]@Polygons[[1]]@coords[,1],sp@polygons[[1]]@Polygons[[1]]@coords[,2])
if(hole && (length(sp@polygons[[1]]@Polygons)>1))
{
#pts in holes
for(k in 2:length(sp@polygons[[1]]@Polygons))
{
logicalPoly=logicalPoly-point.in.polygon(PointsSpatiaux$x,PointsSpatiaux$y,sp@polygons[[1]]@Polygons[[k]]@coords[,1],sp@polygons[[1]]@Polygons[[k]]@coords[,2])
}
}
return(pts[logicalPoly!=0,])
}
##################################################################
#' searchNODcrit
#'
#' @details description, a paragraph
#' @param qProb probability vector used to generate quantile values
#' @param le level index
#' @param zk list of zonings
#' @param criterion list of criteria
#' @param cost list of costs
#' @param costL list of per label costs
#' @param nz list of numbers of zones
#'
#' @return a list with components
#'\describe{
#' \item{ind}{index of last level zoning that has the higher criterion value}
#' \item{critList}{criterion value corresponding to best last level zoning}
#' \item{costlist}{cost value corresponding to best last level zoning}
#' \item{costLlist}{cost per label value corresponding to best last level zoning}
#' \item{nzList}{number of zones of best last level zoning}
#' \item{nq}{lenght of quantile vector}
#' }
#'
#' @keywords internal
#'
#' @examples
#' qProb=c(0.1,0.2);criti=correctionTree(qProb,mapTest)
#' res=geozoning:::searchNODcrit(qProb,2,criti$zk,criti$criterion,criti$cost,criti$costL,criti$nz)
searchNODcrit=function(qProb,le,zk,criterion,cost,costL,nz)
##################################################################
{
# zk: list of zonings
# le: list index
# crit: list of criteria
critf=unlist(criterion[[le]])
costf=unlist(cost[[le]])
costfL=unlist(costL[[le]])
nzf=unlist(nz[[le]])
# check for degenerated zonings
# number of zone labels < number of quantiles+1
# if labels start at 1
nq0=length(qProb)+1
# compute number of labs for each solution
nqf=c()
ind=list()
bestcrit=list()
bestcost=list()
bestcostL=list()
bestnz=list()
nq=list()
lk=length(zk[[le]])
kk=1:lk
for (ilab in 1:lk){
u=unique(zk[[le]][[ilab]]$lab)
lu=length(u)
nqf=c(nqf,lu)
}
while (nq0>1){
maskNOD=(nqf==nq0)
critq=critf[maskNOD]
if(!is.na(critq) && (length(critq)>0)){
ii=which(critq == max(critq)) #search from the best non degenerated ones
critmax=critq[ii]
# find original index in critf vector
k0=kk[maskNOD]
jj=k0[ii]
labq=paste("q",nq0-1,sep="")
ind[[labq]]=jj
bestcrit[[labq]]=critf[jj]
bestcost[[labq]]=costf[jj]
bestcostL[[labq]]=costfL[jj]
bestnz[[labq]]=nzf[jj]
nq[[labq]]=nq0-1
}
nq0=nq0-1
}#end while
return(list(ind=ind,critList=bestcrit,costList=bestcost,costLList=bestcostL,nzList=bestnz,nq=nq))
}
##################################################################
#' searchNODcrit1
#'
#' @details description, a paragraph
#' @param qProb probability vector used to generate quantile values
#' @param crit result of call to correctionTree (with SAVE=TRUE)
#'
#' @return a list with components
#'\describe{
#' \item{ind}{index of last level zoning that has the higher criterion value}
#' \item{critList}{criterion value corresponding to best last level zoning}
#' \item{costlist}{cost value corresponding to best last level zoning}
#' \item{costLlist}{cost per label value corresponding to best last level zoning}
#' \item{nzList}{number of zones of best last level zoning}
#' }
#' @export
#'
#' @examples
#' data(mapTest)
#' qProb=c(0.1,0.2,0.4);criti=correctionTree(qProb,mapTest) # 2 zonings at last level
#' res=searchNODcrit1(qProb,criti)# best one is frist element of last level
searchNODcrit1=function(qProb,crit)
##################################################################
{
# zk: list of zonings
# le: list index
zk=crit$zk
cr=crit$criterion
cost=crit$cost
costL=crit$costL
nz=crit$nz
le=length(zk)
# crit: list of criteria
critf=unlist(cr[[le]])
costf=unlist(cost[[le]])
costfL=unlist(costL[[le]])
nzf=unlist(nz[[le]])
# check for degenerated zonings
# number of zone labels < number of quantiles+1
# if labels start at 1
nq0=length(qProb)+1
# compute number of labs for each solution
nqf=c()
ind=list()
bestcrit=list()
bestcost=list()
bestcostL=list()
bestnz=list()
nq=list()
lk=length(zk[[le]])
kk=1:lk
for (ilab in 1:lk){
u=unique(zk[[le]][[ilab]]$lab)
lu=length(u)
nqf=c(nqf,lu)
}
while (nq0>1){
maskNOD=(nqf==nq0)
critq=critf[maskNOD]
if(!is.na(critq) && (length(critq)>0)){
ii=which(critq == max(critq)) #search from the best non degenerated ones
critmax=critq[ii]
# find original index in critf vector
k0=kk[maskNOD]
jj=k0[ii]
labq=paste("q",nq0-1,sep="")
ind[[labq]]=jj
bestcrit[[labq]]=critf[jj]
bestcost[[labq]]=costf[jj]
bestcostL[[labq]]=costfL[jj]
bestnz[[labq]]=nzf[jj]
nq[[labq]]=nq0-1
}
nq0=nq0-1
}#end while
return(list(ind=ind,critList=bestcrit,costList=bestcost,costLList=bestcostL,nzList=bestnz,nq=nq))
}
##################################################################
#' normDistMat
#'
#' @details normalize distance matrix so that diagonal is equal to 1
#' @param matDistanceCorr corrected distance matrix using pErr, result of calDistance
#' @param optiCrit criterion choice
#'
#' @return a normalized distance matrix
#'
#' @export
#'
#' @examples
#' # load test map with simulated data
#' data(mapTest)
#' # load zoning results from test file
#' data(resZTest)
#' K=resZTest
#' resD = calDistance(typedist=1,mapTest$krigData,K$listZonePoint,K$zoneN,
#' mapTest$krigSurfVoronoi,K$meanZone,pErr=0.9)
#' normDistMat(resD$matDistanceCorr,2)
normDistMat=function(matDistanceCorr,optiCrit)
##################################################################
{
# other values of optiCrit not managed
#
if(optiCrit==4||optiCrit==6)
normMD=distanceNormalisationSqrt(matDistanceCorr)
if(optiCrit==2)
normMD=distanceNormalisationSum(matDistanceCorr)
return(normMD)
}
##################################################################
#' linesC
#'
#' @details add contour Lines to plot
#' @param listContour list of contour lines
#' @param col line color
#'
#' @return an empty value
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' plot(mapTest$boundary)
#' geozoning:::linesC(cL,col="black")
linesC = function(listContour,col="blue")
##################################################################
{
for (i in 1:length(listContour)){
lines(listContour[[i]],col=col)
}
return()
}
##################################################################
#' interCB
#'
#' @details generates SpatialPolygons object corresponding to intersection of contour with boundary, must be within SpatialPolygons given in envel argument
#' @param co contour line
#' @param step map grid resolution
#' @param bd map boundary
#' @param envel envelope
#' @param disp info level (0-no info, 1- add lines to plot)
#' @importFrom sp Polygon
#' @importFrom sp plot
#' @return a SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' pG=polyToSp2(sp::Polygon(mapTest$boundary)) #SpatialPolygons corresponding to map boundary
#' cL=contourAuto(list(),mapTest$step,mapTest$xsize,mapTest$ysize,
#' mapTest$krigGrid,c(5,7),mapTest$boundary)
#' ps = geozoning:::interCB(cL[[8]],mapTest$step,mapTest$boundary,pG)#envelope is the whole map
#' sp::plot(pG)
#' sp::plot(ps,col="red",add=TRUE)
interCB = function(co,step,bd=list(x=c(0,0,1,1,0),y=c(0,1,1,0,0)),envel,disp=0)
##############################################################################
{
#returns spatial polygon corresponding to intersection of contour with boundary
#
polygoneGlobal=SpatialPolygons(list(Polygons(list(Polygon(bd, hole = FALSE)), "1")))
contourL = ContourLines2SLDF(list(co))
polyBuff=gBuffer(contourL,width=0.0001*step)
polyDiff=gDifference(polygoneGlobal,polyBuff)
recupPoly=separationPoly(polyDiff)
ler=length(recupPoly)
if(ler<2) return(NULL) # intersection=boundary -> degenerate contour
sp1=recupPoly[[1]] #
sp2=recupPoly[[2]] #
# keep the smallest one that is within the envelope
sp=sp2
if (gContains(envel,sp1))
{
sp=sp1
if (gContains(envel,sp2) & (gArea(sp2)<gArea(sp1))) sp=sp2
}
if(disp) linesSp(sp)
return(sp)
}
##################################################################
#' getNq
#'
#' @details determine size of quantile in result from correctionTree
#' @param critList component critList of result from correctionTree
#'
#' @return a vector with the size of quantile vectors for each zoning corresponding to critList
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7,
#' # saving only best last level results
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=FALSE)
#' geozoning:::getNq(criti$critList)
getNq=function(critList)
##################################################################
{
n=names(critList)
for (qq in 1:length(critList))
{
nq=sapply(strsplit(n[qq],"q"),function(x){return(x[2])})
nq=as.numeric(nq)
print(paste("criterion=",round(critList[[qq]][1],2),"nq=",nq))
}
return(nq)
}
##################################################################
#' addContour
#'
#' @details add contour lines to plot
#' @param map object returned by function genMap
#' @param val quantile value vector
#' @param col color parameter
#' @param super if TRUE add to existing plot lines coresponding to contour, if FALSE plot boundary and add lines
#'
#' @return void
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' addContour(mapTest,c(5,7),super=FALSE)
addContour=function(map,val,col="blue",super=TRUE)
##################################################################
{
for ( v in val)
{
cL=list()
cL = contourAuto(cL,map$step,map$xsize,map$ysize,map$krigGrid,v,map$boundary)
lc = length(cL)
if (lc >0)
{
sp=list()
if (!super) plot(map$boundary,type="l") # new plot
for (k in 1:lc)
lines(cL[[k]],col=col) # add contour to existing plot
}
}
return()
}
##########################################################################
#' extractionPoly
#'
#' @details extract all elements from SpatialPolygons, holes and full polygons are handled equally
#' @param polyTot SpatialPolygons
#'
#' @return a list of SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' ZK=initialZoning(qProb=c(0.2,0.4,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' geozoning:::extractionPoly(Z[[5]]) # returns 2 SpatialPolygons
extractionPoly=function(polyTot)
###########################################################################
{
inPoly=polyTot@polygons[[1]]@Polygons
listPolyExtract=list()
for(i in (1:length(inPoly)))
{
listPolyExtract[[i]]=SpatialPolygons(list(Polygons(list(inPoly[[i]]), "1")))
}
return(listPolyExtract)
}
###################################
#' plotVario
#'
#' @details plot empirical variogram for model and data in map (raw data plus kriged data)
#' @param map object returned by function genMap
#' @param ylim range of y axis
#'
#' @return a plot
#'
#' @export
#' @importFrom RandomFields RFempiricalvariogram
#'
#' @examples
#' data(mapTest)
#' plotVario(mapTest)
plotVario=function(map,ylim=NULL)
###################################
{
modelGen=map$modelGen
data=map$rawData #raw data
dataK=map$krigData #kriged data
empvario=RFempiricalvariogram(data=data)
empvarioK=RFempiricalvariogram(data=dataK)
if(!is.null(modelGen))
{
plot(empvario,model=modelGen,ylim=ylim)
#kriged variogram
plot(empvarioK,model=modelGen,ylim=ylim,col="blue")
}
else plot(empvario)
}
#########################
#' costLab
#'
#' @details description, a paragraph
#' @param K zoning object, as returned by the calNei function
#' @param map object returned by genMap function
#'
#' @return the sum of per label costs
#'
#' @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)
#' K=criti$zk[[1]][[1]] # initial zoning
#' costLab(K,mapTest) #identical to criti$costL[[1]][[1]]
costLab=function(K,map)
#########################
{
#number of labels
uni=unique(K$lab)
nL=length(uni)
nZ=length(K$lab)
listZonePoint=K$listZonePoint
tabVal=map$krigData
surfVoronoi=map$krigSurfVoronoi
#find which zones are assigned to each label
zlab=list()
vZ=1:nZ
for (ilab in uni){
mask=which(K$lab==ilab)
zlab[[ilab]]=vZ[mask]
}
# compute costL (per label)
res = SigmaL2(zlab,listZonePoint,tabVal,surfVoronoi)
return(res$cL)
}
################################################################
#' SigmaL2
#'
#' @details compute overall mean and variance of all zones for each label plus sum of them for all labels
#' @param zlab list with zone numbers for each zone label
#' @param listZonePoint list of indices of data points within zones, result of call to \code{\link{calNei}}
#' @param tabVal SpatialPointsDataFrame containing data values
#' @param surfVoronoi Surfaces of the Voronoi polygons corresponding to data pts
#'
#' @return a list with components
#' \describe{
#' \item{cL}{weighted (with Voronoi surfaces) average of per label variances}
#' \item{SigmaL2}{vector of per label variances}
#' \item{SL}{vector of per label Voronoi surfaces}
#' \item{mL}{vector of weighted (with Voronoi surfaces) per label average values}
#' \item{voroLab}{vector of per label data}
#' }
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7
#' # save initial zoning and last level zonings
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=TRUE)
#' K=criti$zk[[2]][[1]]
#' uni=unique(K$lab)
#' zlab=sapply(uni,function(x){(1:length(K$lab))[K$lab==x]})
#' sig=SigmaL2(zlab,K$listZonePoint,mapTest$krigData,mapTest$krigSurfVoronoi)
#' }
SigmaL2=function(zlab,listZonePoint,tabVal,surfVoronoi)
################################################################
{
# zlab: list with zone numbers for each zone label
# compute overall mean and std of all zones for each label
# and overall cost per label
# remove empty labels from zlab
zlabNULL=sapply(zlab,is.null)
zlab=zlab[!zlabNULL]
nL=length(zlab)
mL=c()
SL=c()
vLab=list()
voroLab=list()
#first compute mean
for (k in 1:nL){
zlabK=zlab[[k]]
vLabK=c()
voroLabK=c()
for (j in zlabK){
vLabK=c(vLabK,tabVal@data[listZonePoint[[j]],1]) #all data values for zones with label k
voroLabK=c(voroLabK,surfVoronoi[listZonePoint[[j]]]) # corresp. voronoi surfaces
}
vLab[[k]]=vLabK
voroLab[[k]]=voroLabK
SL[k]= sum(voroLabK)
mL[k]=sum(vLabK*voroLabK)/SL[k]
}
# then compute sd
sigmaL2=rep(0,nL)
for (k in 1:nL){
sigmaL2[k]=sum(((vLab[[k]]-mL[k])^2)*voroLab[[k]])/SL[k]
}
cL=sum(sigmaL2*SL)/sum(SL)
return(list(cL=cL,sigmaL2=sigmaL2,SL=SL,mL=mL,voroLab=voroLab))
}
################################################################
#' meanL
#'
#' @details compute overall mean of all zones for each label
#' @param zlab list with zone numbers for each zone label
#' @param listZonePoint list of indices of data points within zones, result of call to \code{\link{calNei}}
#' @param tabVal SpatialPointsDataFrame containing data values
#' @param surfVoronoi Surfaces of the Voronoi polygons corresponding to data pts
#'
#' @return a list with components
#' \describe{
#' \item{mL}{vector of weighted (with Voronoi surfaces) per label average values}
#' \item{SL}{vector of per label Voronoi surfaces}
#' }
#'
#' @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)
#' K=criti$zk[[2]][[1]]
#' uni=unique(K$lab)
#' zlab=sapply(uni,function(x){(1:length(K$lab))[K$lab==x]})
#' resL=meanL(zlab,K$listZonePoint,mapTest$krigData,mapTest$krigSurfVoronoi)
meanL=function(zlab,listZonePoint,tabVal,surfVoronoi)
################################################################
{
# zlab: list with zone numbers for each zone label
nL=length(zlab)
mL=c()
SL=c()
for (k in 1:nL){
zlabK=zlab[[k]]
vLabK=c()
voroLabK=c()
for (j in zlabK){
vLabK=c(vLabK,tabVal@data[listZonePoint[[j]],1]) #all data values for zones with label k
voroLabK=c(voroLabK,surfVoronoi[listZonePoint[[j]]]) # corresp. voronoi surfaces
}
SL[k]= sum(voroLabK)
mL[k]=sum(vLabK*voroLabK)/SL[k]
}
return(list(mL=mL,SL=SL))
}
######################################
#' meansdSimu
#'
#' @details computes mean and standard deviation of a set of map simulations
#' @param vseed list of simulation seeds
#' @param krig type of kriging (1-variogram model-based, 2-inverse distance-based)
#'
#' @return a matrix with as many rows as simulations, and 4 columns, the first two columns give mean and standard deviation of generated raw data, the last two columns give mean and standard deviation of kriged data
#'
#' @export
#'
#' @examples
#' \donttest{
#' meansdSimu(c(1,2))
#' }
meansdSimu=function(vseed=NULL,krig=2)
######################################
{
m=c()
km=m
sdd=m
sdkd=m
for (seed in vseed){
map=genMap(DataObj=NULL,seed=seed,krig=krig,disp=0)
v=map$rawData
v=v@data[,1]
kv=map$krigData
kv=kv@data[,1]
m=c(m,mean(v))
sdd=c(sdd,sd(v))
km=c(km,mean(kv))
sdkd=c(sdkd,sd(kv))
}
mat=cbind(m,km,sdd,sdkd)
rownames(mat)=paste(vseed)
return(mat)
}
######################################
#' meanvarsimu
#'
#' @details computes mean and standard deviation of a set of map simulations
#' @param map object generated by genMap
#'
#' @return a vector with 4 elements, the first two give mean and standard deviation of generated raw data, the last two give mean and standard deviation of kriged data
#'
#' @export
#'
#' @examples
#' \donttest{
#' meanvarSimu(mapTest)
#' }
meanvarSimu=function(map)
######################################
{
v=map$rawData
v=v@data[,1]
kv=map$krigData
kv=kv@data[,1]
m=mean(v)
km=mean(kv)
sdd=sd(v)
sdkd=sd(kv)
vec=c(m,km,sdd,sdkd)
names(vec)=c("raw mean","kriged mean","raw sd","kriged sd")
return(vec)
}
###################
#' valZ
#'
#' @details sorts zones according to attribute mean value
#' @param map map object returned by genMap function
#' @param K zoning object (such as returned by calNei function)
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a list with components
#' \describe{
#' \item{val}{list with vector of data values for each zone, zones are sorted by increasing mean values}
#' \item{ord}{order of zones sorted by increasing mean values}
#' }
#'
#' @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)
#' K=criti$zk[[2]][[1]]
#' valZ(mapTest,K)
valZ=function(map,K)
###################
{
tab=map$krigData
pt=K$listZonePoint
mu=K$meanZone
ord=order(mu)
val=list()
k=0
for(ii in ord){
k=k+1
val[[k]]=tab[[1]][pt[[ii]]]
}
return(list(val=val,ord=ord))
}
###################
#' orderZ
#'
#' @details sorts zones according to ord vector
#' @param Z zoning geometry
#' @param ord sorting order
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a zoning geometry (list of SpatialPolygons)
#'
#' @export
#'
#' @examples
#' \donttest{
#' map=genMap(DataObj=NULL,seed=40,disp=FALSE,krig=1,Vnugget=1.2,typeMod="Gau")
#' qProb=c(0.275,0.8)
#' criti=correctionTree(qProb,map,LASTPASS=FALSE)
#' res=searchNODcrit1(qProb,criti)
#' b=res$ind[[1]][1]
#' K=criti$zk[[2]][[b]]
#' Z=K$zonePolygone
#' plotZ(Z)
#' ord=valZ(map,K)$ord
#' Z=orderZ(Z,ord)
#' plotZ(Z)
#' }
###################
orderZ=function(Z,ord)
{
k=0
for (iZ in ord){
k=k+1
Z=setId(Z,iZ,k)
}
return(Z)
}
#############################
#' superLines
#'
#' @details converts boundary (list of x and y pts) into Spatial Lines
#' @param boundary list, contains x and y coordinates of map boundaries
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a SpatialLines object
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' superL=superLines(mapTest$boundary)
#' sp::plot(superL)
superLines=function(boundary)
#############################
{
#transform boundary into SpatialLines
boundary = data.frame(boundary)
sp::coordinates(boundary)=~x+y
bl=Line(coordinates(boundary))
bspl=SpatialLines(list(Lines(list(bl),'1')))
bdLines = bspl@lines[[1]]@Lines[[1]]
listBdLines=list(Lines(list(Line(bdLines@coords[1:2,])),'1'))
for (i in 2:(length(bdLines@coords)/2-1))
{
listBdLines[[i]] = Lines(list(Line(bdLines@coords[i:(i+1),])),paste(i))
}
SuperLines = SpatialLines(listBdLines)
return(SuperLines)
}
####################
#' r2
#'
#' @details adjusted R2
#' @param reslm result of a call to lm
#'
#' @return the adjusted r-square of the lm model
#'
#' @export
#' @importFrom stats anova dist lm quantile sd
#'
#' @examples
#' # not run
r2=function(reslm)
####################
{
s2T <- sum(anova(reslm)[[2]]) / sum(anova(reslm)[[1]])
MSE <- anova(reslm)[[3]][2]
adj.R2 <- (s2T - MSE) / s2T
return(adj.R2)
}
hazaeljones/geozoning documentation built on May 30, 2019, 3:06 p.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.
###################################################################
#' getCoords
#'
#' @details get SpatialPolygons coordinates
#' @param sp SpatialPolygons
#' @param k polygon number
#'
#' @return some coordinates
#'
#' @keywords internal
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' geozoning:::getCoords(Z[[1]])
getCoords=function(sp,k=1)
##################################################################
{
coords=sp@polygons[[1]]@Polygons[[k]]@coords
return(coords)
}
##################################################################
#' spToSL
#'
#' @details tranform SpatialPolygons into SpatialLines
#' @param sp SpatialPolygons
#'
#' @return a SpatialLines
#'
#' @export
#' @importFrom sp SpatialLines Lines
#'
#' @examples
#' data(resZTest)
#' K=resZTest
#' Z=K$zonePolygone
#' spToSL(Z[[5]])
spToSL=function(sp)
##################################################################
{
co = getCoords(sp)
li = Lines(list(Line(co)),ID="1")
lis = SpatialLines(list(li))
return(lis)
}
##################################################################
#' contourArea
#'
#' @details area corresponding to closed contour line
#' @param co contour line
#'
#' @return the area within the contour line
#'
#' @export
#' @importFrom sp SpatialPolygons SpatialPointsDataFrame Polygons Polygon
#' @importFrom maptools ContourLines2SLDF
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' contourArea(cL[[8]])
contourArea=function(co)
##################################################################
{
contour = ContourLines2SLDF(list(co))
coords = contour@lines[[1]]@Lines[[1]]@coords
poly=Polygon(coords)
polys=Polygons(list(poly),"id")
contourSp = SpatialPolygons(list(polys))
surface = gArea(contourSp)
return(surface)
}
##################################################################
#' listContourArea
#'
#' @details area of all contour lines in list
#' @param cL list of contour lines
#'
#' @return a list of areas
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' geozoning:::listContourArea(cL)
listContourArea=function(cL)
##################################################################
{
return(sapply(cL,contourArea))
}
##################################################################
#' contourToSpp
#'
#' @details transform contour line into SpatialPolygons
#' @param co contour line (list with contour level and x,y coordinates
#' @param step grid resolution
#'
#' @return a list with components
#'\describe{
#'\item{sp}{SpatialPolygons corresponding to contour line}
#'\item{contour}{SpatialLines corresponding to contour line}
#'\item{polyBuff}{SpatialPolygons corresponding to buffer around contour line}
#'\item{surface}{SpatialPolygons area}
#'}
#'
#' @export
#' @importFrom rgeos gBuffer gArea gCentroid gContains gConvexHull gDifference gDistance gIntersects gWithin
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' contourToSpp(cL[[8]],0.1)
contourToSpp=function(co,step)
##################################################################
{
contour = ContourLines2SLDF(list(co))
coords = contour@lines[[1]]@Lines[[1]]@coords
# attention gBuffer renvoie 2 polygones
polyBuff = gBuffer(contour,width=0.0001*(1/step))
poly=Polygon(coords)
polys=Polygons(list(poly),"id")
contourSp = SpatialPolygons(list(polys))
surface = gArea(contourSp)
return(list(sp=contourSp,contour=contour,polyBuff=polyBuff,surface=surface))
}
##################################################################
#' nPolySp
#' @details number of polygons in SpatialPolygons
#' @param sp SpatialPolygons
#' @return the number of polygons within the current zone
#'
#' @keywords internal
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' print(paste(geozoning:::nPolySp(Z[[2]]),"polygons"))
nPolySp =function(sp)
##################################################################
{
return(length(sp@polygons[[1]]@Polygons))
}
##################################################################
#' holeSp
#'
#' @details number of holes in SpatialPolygons
#' @param sp SpatialPolygons
#'
#' @return the number of holes within sp
#'
#' @export
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' holeSp(Z[[5]]) #zone 5 has 1 hole
holeSp = function(sp)
##################################################################
{
le=nPolySp(sp)
nh=0
if (le == 0) return(0)
for (i in 1:le)
{
poly=sp@polygons[[1]]@Polygons[[i]]
if(poly@hole) nh=nh+1
}
return (nh)
}
##################################################################
#' maxDistSP
#'
#' @details maximum distance within kth polygon of current zone
#' @param sp SpatialPolygons
#' @return the maximum distance within sp
#'
#' @keywords internal
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' geozoning:::maxDistSP(Z[[5]])
maxDistSP=function(sp)
##################################################################
{
return(max(dist(sp@polygons[[1]]@Polygons[[1]]@coords)))
}
##################################################################
#' getPolySp
#'
#' @details get the kth polygon of the current SpatialPolygons
#' @param sp SpatialPolygons object
#' @param k polygon number
#'
#' @return a polygon (object of class Polygon)
#'
#' @keywords internal
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' sp=Z[[5]]
#' P1=geozoning:::getPolySp(sp,1)
#' P2=geozoning:::getPolySp(sp,2) # second polygon is a hole
#' plot(P1@coords,type="l")
#' lines(P2@coords,type="l",col="blue")
getPolySp = function(sp,k=1)
##################################################################
{
p=sp@polygons[[1]]@Polygons[[k]]
return(p)
}
##################################################################
#' polyToSp2
#'
#' @details transforms polygon into SpatialPolygons
#' @param p polygon
#'
#' @return a SpatialPolygons
#'
#' @export
#'
#' @examples
#' ZK=initialZoning(qProb=c(0.4,0.2,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' sp=Z[[5]]
#' P1=geozoning:::getPolySp(sp,1)
#' sph=polyToSp2(P1)
#' plotZ(Z)
#' sp::plot(sph,col="blue",lwd=2,add=TRUE)
polyToSp2=function(p)
##################################################################
{
sp=SpatialPolygons(list(Polygons(list(p),1:1)))
return(sp)
}
##################################################################
#' lineToSp
#'
#' @details transform closed line into SpatialPolygons
#' @param lin list with x and y line coordinates
#'
#' @return a SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' lin=data.frame(x=cL[[8]]$x,y=cL[[8]]$y)
#' sp=geozoning:::lineToSp(lin)
lineToSp=function(lin)
##################################################################
{
sp=SpatialPolygons(list(Polygons(list(Polygon(lin,hole = FALSE)), "1")))
return(sp)
}
##################################################################
#' spnorm
#'
#' @details normalize Polygon according to border limits
#' @param sp object of class Polygons
#' @param boundary list with x and y components, used to normalize sp
#'
#' @return a list with components
#' \describe{
#' \item{pn}{normalized Polygon}
#' \item{boundaryn}{normalized boundary}
#' }
#'
#' @keywords internal
#'
#' @examples
#' # shape1: result of call to readS on shapefile
#' z=geozoning::shape1
#' bb=list(x=z@bbox[1,],y=z@bbox[2,])
#' p=z@polygons
#' p1=p[[1]]
#' P1=p1@Polygons[[1]]
#' NP1=geozoning:::spnorm(P1,bb)$pn
#' Nbb=geozoning:::spnorm(P1,bb)$boundaryn
#' plot(NP1@coords,xlim=Nbb$x,ylim=Nbb$y)
spnorm = function(sp, boundary)
##################################################################
{
xmin=min(boundary$x)
xmax=max(boundary$x)
ymin=min(boundary$y)
ymax=max(boundary$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
if (abs(ymax-ymin) < 1e-6) return(NULL)
x= sp@coords[,1]
y = sp@coords[,2]
x = (x-xmin)/(xmax-xmin)
y = (y-ymin)/(ymax-ymin)
pn = Polygon(cbind(x,y))
bx=boundary$x
by=boundary$y
bx = (bx-xmin)/(xmax-xmin)
by = (by-ymin)/(ymax-ymin)
bn = list(x=bx,y=by)
return(list(pn=pn, boundaryn=bn))
}
##################################################################
#' normalize data coordinates and border
#'
#' @details normalize boundary between 0 and 1 and data coordinates accordingly
#' @param data data frame with x and y components
#' @param bd boundary (list with x and y components)
#'
#' @return a list with components
#'\describe{
#' \item{dataN}{normalized data}
#' \item{boundaryN}{normalized boundary}
#' \item{xmin}{minimum vaue of x within boundary}
#' \item{xmax}{maximum vaue of x within boundary}
#' \item{ymin}{minimum vaue of y within boundary}
#' \item{ymax}{maximum vaue of y within boundary}
#' }
#'
#' @keywords internal
#'
#' @examples
#' x=runif(100, min=0, max=1)
#' y=runif(100, min=0.2, max=1.7)
#' range(x) # not [0,1]
#' tabData=data.frame(x=x,y=y)
#' bd=list(x=c(0,0,1,1,0), y=c(0.2,1.7,1.7,0.2,0.2))
#' res=geozoning:::datanorm(tabData,bd)
#' apply(res$dataN,2,range)#
datanorm = function(data, bd)
##################################################################
{
xmin=min(bd$x)
xmax=max(bd$x)
ymin=min(bd$y)
ymax=max(bd$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
if (abs(ymax-ymin) < 1e-6) return(NULL)
x = data$x
y = data$y
x = (x-xmin)/(xmax-xmin)
y = (y-ymin)/(ymax-ymin)
data$x = x
data$y = y
#normalize border
bx=bd$x
by=bd$y
bx = (bx-xmin)/(xmax-xmin)
by = (by-ymin)/(ymax-ymin)
bdn = list(x=bx,y=by)
return(list(dataN=data, boundaryN=bdn,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax))
}
##################################################################
#' normalize data coords
#'
#' @details normalize data coordinates between 0 and 1 with different ratios for x and y
#' @param data frame with x and y components
#'
#' @return a normalized data frame
#'
#' @keywords internal
#'
#' @examples
#' nPoints=500
#' x=runif(nPoints, min=0, max=1)
#' y=runif(nPoints, min=0, max=1)
#' range(x) # not [0,1]
#' tabData=data.frame(x=x,y=y)
#' tabData=geozoning:::datanormXY(tabData) # x,y ranges are now [0,1]
datanormXY = function(data)
##################################################################
{
xmin=min(data$x)
xmax=max(data$x)
ymin=min(data$y)
ymax=max(data$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
if (abs(ymax-ymin) < 1e-6) return(NULL)
x = data$x
y = data$y
x = (x-xmin)/(xmax-xmin)
y = (y-ymin)/(ymax-ymin)
data$x = x
data$y = y
return(data)
}
##################################################################
#' normalize data coords with same ratio (for non square field)
#'
#' @details normalize x between 0 and 1, y and boundary with same ratio
#' @param data frame with x and y components
#' @param bd list with x and y components
#'
#' @return a list with components
#'\describe{
#' \item{dataN}{normalized data}
#' \item{boundaryN}{normalized boundary}
#' \item{ratio}{normalizing ratio}
#' \item{xmin}{minimum value of x within boundary}
#' \item{xmax}{maximum value of x within boundary}
#' \item{ymin}{minimum value of y within boundary}
#' \item{ymax}{maximum value of y within boundary}
#' }
#'
#' @export
#'
#' @examples
#' x=runif(100, min=0, max=1)
#' y=runif(100, min=0.2, max=1.7)
#' range(x) # not [0,1]
#' tabData=data.frame(x=x,y=y)
#' bd=list(x=c(0,0,1,1,0), y=c(0.2,1.7,1.7,0.2,0.2))
#' res=datanormX(tabData,bd)
#' apply(res$dataN,2,range)# x range is now [0,1], not y range
#' res$ratio # normalization ratio
datanormX = function(data,bd)
##################################################################
{
xmin=min(data$x)
xmax=max(data$x)
ymin=min(data$y)
ymax=max(data$y)
if (abs(xmax-xmin) < 1e-6) return(NULL)
x = data$x
y = data$y
ratio=xmax-xmin
x = (x-xmin)/ratio
y = (y-ymin)/ratio
data$x = x
data$y = y
#normalize boarder
bx=bd$x
by=bd$y
bx = (bx-xmin)/ratio
by = (by-ymin)/ratio
bdn = list(x=bx,y=by)
return(list(dataN=data, boundaryN=bdn,ratio=ratio,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax))
}
##################################################################
#' normSize
#'
#' @details normalize thresholds for small zone detection and no grow zone, considering mapo boundary
#' @param boundaryN normalized map boundary
#' @param minSize minimum size threshold
#' @param minSizeNG no grow size threshold
#'
#' @return a list with components
#' \describe{
#' \item{minSize}{normalized minimum size threshold}
#' \item{minSizeNG}{normalized no grow size threshold}
#' }
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' resT=normSize(mapTest$boundary,0.012,0.001)#normalize thresholds relatively to map boundary area
normSize=function(boundaryN,minSize,minSizeNG)
##################################################################
{
# normalize threshold for small zone detection
coords=boundaryN
poly=Polygon(coords)
polys=Polygons(list(poly),"id")
boundarySp = SpatialPolygons(list(polys))
boundaryArea=gArea(boundarySp)
minSize=minSize/boundaryArea
minSizeNG=minSizeNG/boundaryArea
print(paste("after standardization minSize=",minSize,"minSizeNG=",minSizeNG))
return(list(minSize=minSize,minSizeNG=minSizeNG))
}
##################################################################
#' genQseq
#'
#' @details description, a paragraph
#' @param qProb probability vector used to generate quantile values
#' @param K zoning object, as returned by the calNei function
#' @param map object returned by function genMap
#' @param i1 current zone index
#' @param i2 englobing zone index
#' @param LEQ length of quantile sequence
#' @param MAXP maximum shift from center for quantile sequence
#' @param disp 0: no info, 1: some info
#'
#' @return a plot
#'
#' @export
#' @examples
#' data(mapTest)
#' qProb=c(0.4,0.7)
#' ZK=initialZoning(qProb,mapTest)
#' K=ZK$resZ
#' print(K$lab)
#' genQseq(qProb,K,mapTest,1,2) # from label 3 to label 2
genQseq = function(qProb,K,map,i1,i2,LEQ=5,MAXP=0.1,disp=0)
##################################################################
{
# i1 = current zone index, i2= englobing zone index
#
lab1 = K$lab[i1]
lab2 = K$lab[i2]
nq = length(qProb)
inc = TRUE # increase quantile value
if (lab1 >= lab2) inc = FALSE # decrease quantile value
valRef= quantile(map$krigGrid,na.rm=TRUE,prob=qProb)
# find quantile value corresponding to zone i1
q1min = max(1,lab1-1)
q1max = min(lab1,nq)
if(inc)
q1=q1max
else
q1=q1min
prob1 = qProb[q1]
q1= quantile(map$krigGrid,na.rm=TRUE,prob=prob1)
# increase or decrease quantile value
if (inc)
{
q2 = quantile(map$krigGrid,na.rm=TRUE,prob= min(prob1 + MAXP, 0.99))
Qseq =seq(q1,q2,length.out=LEQ+1)
}
else
{
q3 = quantile(map$krigGrid,na.rm=TRUE,prob= max(prob1 - MAXP, 0.01))
Qseq = seq(q1,q3,length.out=LEQ+1)
}
Qseq=Qseq[-1]
Qseq=sort(unique(Qseq))
if (disp>0)
{
print("Qseq=")
print(Qseq)
}
return(Qseq)
}
##################################################################
#' checkContour
#'
#' @details check admissibility for contour line: surface >minSizeNG and refPoint close enough
#' @param contourSp SpatialPolygons corresponding to closed contour line
#' @param step grid resolution
#' @param refPoint referene point
#' @param minSizeNG zone area threshold under which a zone is not admissible
#'
#' @return Null if contour is not admissible or a list with components
#' \describe{
#' \item{contourSp}{SpatialPolygons corresponding to admissible contour}
#' \item{}{polyBuff}{SpatialPolygons corresponding to gBuffer around admissible contour}
#' }
#' @importFrom sp Polygon
#' @importFrom rgeos plot
#' @importFrom rgeos gCentroid
#' @export
#'
#' @examples
#' data(mapTest)
#' cL=contourAuto(list(),mapTest$step,mapTest$xsize,mapTest$ysize,
#' mapTest$krigGrid,c(5,7),mapTest$boundary)
#' pG=polyToSp2(sp::Polygon(mapTest$boundary)) #SpatialPolygons corresponding to map boundary
#' rgeos::plot(pG)
#' sp8 = contourToSpp(cL[[8]],0.1)$sp
#' refPoint = rgeos::gCentroid(sp8)
#' resp=checkContour(sp8,mapTest$step,refPoint)
#' rgeos::plot(resp$contourSp,col="red",add=TRUE)
checkContour = function(contourSp,step,refPoint,minSizeNG=1e-3)
##################################################################
{
# contourSp spatial object
polyBuff = gBuffer(contourSp,width=0.0001*(1/step)) #spatialPolygon
surface = gArea(contourSp)
condi = surface <=minSizeNG || gDistance(gCentroid(polyBuff),refPoint)>=0.1
if (condi)
return (NULL)
else
return (list(contourSp=contourSp,polyBuff=polyBuff))
}
##################################################################
#' cleanSp
#'
#' @details removes from sp polygons that are too small (artefacts of gDifference)
#' @param sp SpatialPolygons
#' @param tol minimum area for removal
#' @return a SpatialPolygons
#'
#' @export
#'
#' @examples
#' # not run
cleanSp = function(sp,tol=1e-5)
##################################################################
{
# number of polygons
n=nPolySp(sp)
polyl=list()
k=0
# remove polygons that are too small (<1e-6)
# (artefacts of gDifference)
for (i in 1:n)
{
poly=sp@polygons[[1]]@Polygons[[i]]
area=poly@area
if (area >= tol)
{
k=k+1
polyl[[k]]=poly
}
}
if (length(polyl)==0) return(NULL)
polys=Polygons(polyl,"id")
spc= SpatialPolygons(list(polys))
return(spc)
}
##################################################################
#' ptInSp
#'
#' @details finds data points in sp
#' @param sp SpatialPolygons
#' @param pts data points
#' @param hole if TRUE also consider points in holes
#'
#' @return a data frame with data points within sp
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' data(resZTest)
#' Z=resZTest$zonePolygone
#' ptsInSp(Z[[5]],mapTest$krigData) # 5 data points within zone 5
ptsInSp=function(sp,pts,hole=FALSE)
##################################################################
{
PointsSpatiaux=SpatialPoints(pts)
logicalPoly=point.in.polygon(PointsSpatiaux$x,PointsSpatiaux$y,sp@polygons[[1]]@Polygons[[1]]@coords[,1],sp@polygons[[1]]@Polygons[[1]]@coords[,2])
if(hole && (length(sp@polygons[[1]]@Polygons)>1))
{
#pts in holes
for(k in 2:length(sp@polygons[[1]]@Polygons))
{
logicalPoly=logicalPoly-point.in.polygon(PointsSpatiaux$x,PointsSpatiaux$y,sp@polygons[[1]]@Polygons[[k]]@coords[,1],sp@polygons[[1]]@Polygons[[k]]@coords[,2])
}
}
return(pts[logicalPoly!=0,])
}
##################################################################
#' searchNODcrit
#'
#' @details description, a paragraph
#' @param qProb probability vector used to generate quantile values
#' @param le level index
#' @param zk list of zonings
#' @param criterion list of criteria
#' @param cost list of costs
#' @param costL list of per label costs
#' @param nz list of numbers of zones
#'
#' @return a list with components
#'\describe{
#' \item{ind}{index of last level zoning that has the higher criterion value}
#' \item{critList}{criterion value corresponding to best last level zoning}
#' \item{costlist}{cost value corresponding to best last level zoning}
#' \item{costLlist}{cost per label value corresponding to best last level zoning}
#' \item{nzList}{number of zones of best last level zoning}
#' \item{nq}{lenght of quantile vector}
#' }
#'
#' @keywords internal
#'
#' @examples
#' qProb=c(0.1,0.2);criti=correctionTree(qProb,mapTest)
#' res=geozoning:::searchNODcrit(qProb,2,criti$zk,criti$criterion,criti$cost,criti$costL,criti$nz)
searchNODcrit=function(qProb,le,zk,criterion,cost,costL,nz)
##################################################################
{
# zk: list of zonings
# le: list index
# crit: list of criteria
critf=unlist(criterion[[le]])
costf=unlist(cost[[le]])
costfL=unlist(costL[[le]])
nzf=unlist(nz[[le]])
# check for degenerated zonings
# number of zone labels < number of quantiles+1
# if labels start at 1
nq0=length(qProb)+1
# compute number of labs for each solution
nqf=c()
ind=list()
bestcrit=list()
bestcost=list()
bestcostL=list()
bestnz=list()
nq=list()
lk=length(zk[[le]])
kk=1:lk
for (ilab in 1:lk){
u=unique(zk[[le]][[ilab]]$lab)
lu=length(u)
nqf=c(nqf,lu)
}
while (nq0>1){
maskNOD=(nqf==nq0)
critq=critf[maskNOD]
if(!is.na(critq) && (length(critq)>0)){
ii=which(critq == max(critq)) #search from the best non degenerated ones
critmax=critq[ii]
# find original index in critf vector
k0=kk[maskNOD]
jj=k0[ii]
labq=paste("q",nq0-1,sep="")
ind[[labq]]=jj
bestcrit[[labq]]=critf[jj]
bestcost[[labq]]=costf[jj]
bestcostL[[labq]]=costfL[jj]
bestnz[[labq]]=nzf[jj]
nq[[labq]]=nq0-1
}
nq0=nq0-1
}#end while
return(list(ind=ind,critList=bestcrit,costList=bestcost,costLList=bestcostL,nzList=bestnz,nq=nq))
}
##################################################################
#' searchNODcrit1
#'
#' @details description, a paragraph
#' @param qProb probability vector used to generate quantile values
#' @param crit result of call to correctionTree (with SAVE=TRUE)
#'
#' @return a list with components
#'\describe{
#' \item{ind}{index of last level zoning that has the higher criterion value}
#' \item{critList}{criterion value corresponding to best last level zoning}
#' \item{costlist}{cost value corresponding to best last level zoning}
#' \item{costLlist}{cost per label value corresponding to best last level zoning}
#' \item{nzList}{number of zones of best last level zoning}
#' }
#' @export
#'
#' @examples
#' data(mapTest)
#' qProb=c(0.1,0.2,0.4);criti=correctionTree(qProb,mapTest) # 2 zonings at last level
#' res=searchNODcrit1(qProb,criti)# best one is frist element of last level
searchNODcrit1=function(qProb,crit)
##################################################################
{
# zk: list of zonings
# le: list index
zk=crit$zk
cr=crit$criterion
cost=crit$cost
costL=crit$costL
nz=crit$nz
le=length(zk)
# crit: list of criteria
critf=unlist(cr[[le]])
costf=unlist(cost[[le]])
costfL=unlist(costL[[le]])
nzf=unlist(nz[[le]])
# check for degenerated zonings
# number of zone labels < number of quantiles+1
# if labels start at 1
nq0=length(qProb)+1
# compute number of labs for each solution
nqf=c()
ind=list()
bestcrit=list()
bestcost=list()
bestcostL=list()
bestnz=list()
nq=list()
lk=length(zk[[le]])
kk=1:lk
for (ilab in 1:lk){
u=unique(zk[[le]][[ilab]]$lab)
lu=length(u)
nqf=c(nqf,lu)
}
while (nq0>1){
maskNOD=(nqf==nq0)
critq=critf[maskNOD]
if(!is.na(critq) && (length(critq)>0)){
ii=which(critq == max(critq)) #search from the best non degenerated ones
critmax=critq[ii]
# find original index in critf vector
k0=kk[maskNOD]
jj=k0[ii]
labq=paste("q",nq0-1,sep="")
ind[[labq]]=jj
bestcrit[[labq]]=critf[jj]
bestcost[[labq]]=costf[jj]
bestcostL[[labq]]=costfL[jj]
bestnz[[labq]]=nzf[jj]
nq[[labq]]=nq0-1
}
nq0=nq0-1
}#end while
return(list(ind=ind,critList=bestcrit,costList=bestcost,costLList=bestcostL,nzList=bestnz,nq=nq))
}
##################################################################
#' normDistMat
#'
#' @details normalize distance matrix so that diagonal is equal to 1
#' @param matDistanceCorr corrected distance matrix using pErr, result of calDistance
#' @param optiCrit criterion choice
#'
#' @return a normalized distance matrix
#'
#' @export
#'
#' @examples
#' # load test map with simulated data
#' data(mapTest)
#' # load zoning results from test file
#' data(resZTest)
#' K=resZTest
#' resD = calDistance(typedist=1,mapTest$krigData,K$listZonePoint,K$zoneN,
#' mapTest$krigSurfVoronoi,K$meanZone,pErr=0.9)
#' normDistMat(resD$matDistanceCorr,2)
normDistMat=function(matDistanceCorr,optiCrit)
##################################################################
{
# other values of optiCrit not managed
#
if(optiCrit==4||optiCrit==6)
normMD=distanceNormalisationSqrt(matDistanceCorr)
if(optiCrit==2)
normMD=distanceNormalisationSum(matDistanceCorr)
return(normMD)
}
##################################################################
#' linesC
#'
#' @details add contour Lines to plot
#' @param listContour list of contour lines
#' @param col line color
#'
#' @return an empty value
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' cL=list()
#' cL=contourAuto(cL,mapTest$step,mapTest$xsize,mapTest$ysize,mapTest$krigGrid,c(5,7),mapTest$boundary)
#' plot(mapTest$boundary)
#' geozoning:::linesC(cL,col="black")
linesC = function(listContour,col="blue")
##################################################################
{
for (i in 1:length(listContour)){
lines(listContour[[i]],col=col)
}
return()
}
##################################################################
#' interCB
#'
#' @details generates SpatialPolygons object corresponding to intersection of contour with boundary, must be within SpatialPolygons given in envel argument
#' @param co contour line
#' @param step map grid resolution
#' @param bd map boundary
#' @param envel envelope
#' @param disp info level (0-no info, 1- add lines to plot)
#' @importFrom sp Polygon
#' @importFrom sp plot
#' @return a SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' pG=polyToSp2(sp::Polygon(mapTest$boundary)) #SpatialPolygons corresponding to map boundary
#' cL=contourAuto(list(),mapTest$step,mapTest$xsize,mapTest$ysize,
#' mapTest$krigGrid,c(5,7),mapTest$boundary)
#' ps = geozoning:::interCB(cL[[8]],mapTest$step,mapTest$boundary,pG)#envelope is the whole map
#' sp::plot(pG)
#' sp::plot(ps,col="red",add=TRUE)
interCB = function(co,step,bd=list(x=c(0,0,1,1,0),y=c(0,1,1,0,0)),envel,disp=0)
##############################################################################
{
#returns spatial polygon corresponding to intersection of contour with boundary
#
polygoneGlobal=SpatialPolygons(list(Polygons(list(Polygon(bd, hole = FALSE)), "1")))
contourL = ContourLines2SLDF(list(co))
polyBuff=gBuffer(contourL,width=0.0001*step)
polyDiff=gDifference(polygoneGlobal,polyBuff)
recupPoly=separationPoly(polyDiff)
ler=length(recupPoly)
if(ler<2) return(NULL) # intersection=boundary -> degenerate contour
sp1=recupPoly[[1]] #
sp2=recupPoly[[2]] #
# keep the smallest one that is within the envelope
sp=sp2
if (gContains(envel,sp1))
{
sp=sp1
if (gContains(envel,sp2) & (gArea(sp2)<gArea(sp1))) sp=sp2
}
if(disp) linesSp(sp)
return(sp)
}
##################################################################
#' getNq
#'
#' @details determine size of quantile in result from correctionTree
#' @param critList component critList of result from correctionTree
#'
#' @return a vector with the size of quantile vectors for each zoning corresponding to critList
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7,
#' # saving only best last level results
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=FALSE)
#' geozoning:::getNq(criti$critList)
getNq=function(critList)
##################################################################
{
n=names(critList)
for (qq in 1:length(critList))
{
nq=sapply(strsplit(n[qq],"q"),function(x){return(x[2])})
nq=as.numeric(nq)
print(paste("criterion=",round(critList[[qq]][1],2),"nq=",nq))
}
return(nq)
}
##################################################################
#' addContour
#'
#' @details add contour lines to plot
#' @param map object returned by function genMap
#' @param val quantile value vector
#' @param col color parameter
#' @param super if TRUE add to existing plot lines coresponding to contour, if FALSE plot boundary and add lines
#'
#' @return void
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' addContour(mapTest,c(5,7),super=FALSE)
addContour=function(map,val,col="blue",super=TRUE)
##################################################################
{
for ( v in val)
{
cL=list()
cL = contourAuto(cL,map$step,map$xsize,map$ysize,map$krigGrid,v,map$boundary)
lc = length(cL)
if (lc >0)
{
sp=list()
if (!super) plot(map$boundary,type="l") # new plot
for (k in 1:lc)
lines(cL[[k]],col=col) # add contour to existing plot
}
}
return()
}
##########################################################################
#' extractionPoly
#'
#' @details extract all elements from SpatialPolygons, holes and full polygons are handled equally
#' @param polyTot SpatialPolygons
#'
#' @return a list of SpatialPolygons
#'
#' @keywords internal
#'
#' @examples
#' data(mapTest)
#' ZK=initialZoning(qProb=c(0.2,0.4,0.7),mapTest)
#' Z=ZK$resZ$zonePolygone
#' geozoning:::extractionPoly(Z[[5]]) # returns 2 SpatialPolygons
extractionPoly=function(polyTot)
###########################################################################
{
inPoly=polyTot@polygons[[1]]@Polygons
listPolyExtract=list()
for(i in (1:length(inPoly)))
{
listPolyExtract[[i]]=SpatialPolygons(list(Polygons(list(inPoly[[i]]), "1")))
}
return(listPolyExtract)
}
###################################
#' plotVario
#'
#' @details plot empirical variogram for model and data in map (raw data plus kriged data)
#' @param map object returned by function genMap
#' @param ylim range of y axis
#'
#' @return a plot
#'
#' @export
#' @importFrom RandomFields RFempiricalvariogram
#'
#' @examples
#' data(mapTest)
#' plotVario(mapTest)
plotVario=function(map,ylim=NULL)
###################################
{
modelGen=map$modelGen
data=map$rawData #raw data
dataK=map$krigData #kriged data
empvario=RFempiricalvariogram(data=data)
empvarioK=RFempiricalvariogram(data=dataK)
if(!is.null(modelGen))
{
plot(empvario,model=modelGen,ylim=ylim)
#kriged variogram
plot(empvarioK,model=modelGen,ylim=ylim,col="blue")
}
else plot(empvario)
}
#########################
#' costLab
#'
#' @details description, a paragraph
#' @param K zoning object, as returned by the calNei function
#' @param map object returned by genMap function
#'
#' @return the sum of per label costs
#'
#' @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)
#' K=criti$zk[[1]][[1]] # initial zoning
#' costLab(K,mapTest) #identical to criti$costL[[1]][[1]]
costLab=function(K,map)
#########################
{
#number of labels
uni=unique(K$lab)
nL=length(uni)
nZ=length(K$lab)
listZonePoint=K$listZonePoint
tabVal=map$krigData
surfVoronoi=map$krigSurfVoronoi
#find which zones are assigned to each label
zlab=list()
vZ=1:nZ
for (ilab in uni){
mask=which(K$lab==ilab)
zlab[[ilab]]=vZ[mask]
}
# compute costL (per label)
res = SigmaL2(zlab,listZonePoint,tabVal,surfVoronoi)
return(res$cL)
}
################################################################
#' SigmaL2
#'
#' @details compute overall mean and variance of all zones for each label plus sum of them for all labels
#' @param zlab list with zone numbers for each zone label
#' @param listZonePoint list of indices of data points within zones, result of call to \code{\link{calNei}}
#' @param tabVal SpatialPointsDataFrame containing data values
#' @param surfVoronoi Surfaces of the Voronoi polygons corresponding to data pts
#'
#' @return a list with components
#' \describe{
#' \item{cL}{weighted (with Voronoi surfaces) average of per label variances}
#' \item{SigmaL2}{vector of per label variances}
#' \item{SL}{vector of per label Voronoi surfaces}
#' \item{mL}{vector of weighted (with Voronoi surfaces) per label average values}
#' \item{voroLab}{vector of per label data}
#' }
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(mapTest)
#' # run zoning with 2 quantiles corresponding to probability values 0.4 and 0.7
#' # save initial zoning and last level zonings
#' criti=correctionTree(c(0.4,0.7),mapTest,SAVE=TRUE)
#' K=criti$zk[[2]][[1]]
#' uni=unique(K$lab)
#' zlab=sapply(uni,function(x){(1:length(K$lab))[K$lab==x]})
#' sig=SigmaL2(zlab,K$listZonePoint,mapTest$krigData,mapTest$krigSurfVoronoi)
#' }
SigmaL2=function(zlab,listZonePoint,tabVal,surfVoronoi)
################################################################
{
# zlab: list with zone numbers for each zone label
# compute overall mean and std of all zones for each label
# and overall cost per label
# remove empty labels from zlab
zlabNULL=sapply(zlab,is.null)
zlab=zlab[!zlabNULL]
nL=length(zlab)
mL=c()
SL=c()
vLab=list()
voroLab=list()
#first compute mean
for (k in 1:nL){
zlabK=zlab[[k]]
vLabK=c()
voroLabK=c()
for (j in zlabK){
vLabK=c(vLabK,tabVal@data[listZonePoint[[j]],1]) #all data values for zones with label k
voroLabK=c(voroLabK,surfVoronoi[listZonePoint[[j]]]) # corresp. voronoi surfaces
}
vLab[[k]]=vLabK
voroLab[[k]]=voroLabK
SL[k]= sum(voroLabK)
mL[k]=sum(vLabK*voroLabK)/SL[k]
}
# then compute sd
sigmaL2=rep(0,nL)
for (k in 1:nL){
sigmaL2[k]=sum(((vLab[[k]]-mL[k])^2)*voroLab[[k]])/SL[k]
}
cL=sum(sigmaL2*SL)/sum(SL)
return(list(cL=cL,sigmaL2=sigmaL2,SL=SL,mL=mL,voroLab=voroLab))
}
################################################################
#' meanL
#'
#' @details compute overall mean of all zones for each label
#' @param zlab list with zone numbers for each zone label
#' @param listZonePoint list of indices of data points within zones, result of call to \code{\link{calNei}}
#' @param tabVal SpatialPointsDataFrame containing data values
#' @param surfVoronoi Surfaces of the Voronoi polygons corresponding to data pts
#'
#' @return a list with components
#' \describe{
#' \item{mL}{vector of weighted (with Voronoi surfaces) per label average values}
#' \item{SL}{vector of per label Voronoi surfaces}
#' }
#'
#' @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)
#' K=criti$zk[[2]][[1]]
#' uni=unique(K$lab)
#' zlab=sapply(uni,function(x){(1:length(K$lab))[K$lab==x]})
#' resL=meanL(zlab,K$listZonePoint,mapTest$krigData,mapTest$krigSurfVoronoi)
meanL=function(zlab,listZonePoint,tabVal,surfVoronoi)
################################################################
{
# zlab: list with zone numbers for each zone label
nL=length(zlab)
mL=c()
SL=c()
for (k in 1:nL){
zlabK=zlab[[k]]
vLabK=c()
voroLabK=c()
for (j in zlabK){
vLabK=c(vLabK,tabVal@data[listZonePoint[[j]],1]) #all data values for zones with label k
voroLabK=c(voroLabK,surfVoronoi[listZonePoint[[j]]]) # corresp. voronoi surfaces
}
SL[k]= sum(voroLabK)
mL[k]=sum(vLabK*voroLabK)/SL[k]
}
return(list(mL=mL,SL=SL))
}
######################################
#' meansdSimu
#'
#' @details computes mean and standard deviation of a set of map simulations
#' @param vseed list of simulation seeds
#' @param krig type of kriging (1-variogram model-based, 2-inverse distance-based)
#'
#' @return a matrix with as many rows as simulations, and 4 columns, the first two columns give mean and standard deviation of generated raw data, the last two columns give mean and standard deviation of kriged data
#'
#' @export
#'
#' @examples
#' \donttest{
#' meansdSimu(c(1,2))
#' }
meansdSimu=function(vseed=NULL,krig=2)
######################################
{
m=c()
km=m
sdd=m
sdkd=m
for (seed in vseed){
map=genMap(DataObj=NULL,seed=seed,krig=krig,disp=0)
v=map$rawData
v=v@data[,1]
kv=map$krigData
kv=kv@data[,1]
m=c(m,mean(v))
sdd=c(sdd,sd(v))
km=c(km,mean(kv))
sdkd=c(sdkd,sd(kv))
}
mat=cbind(m,km,sdd,sdkd)
rownames(mat)=paste(vseed)
return(mat)
}
######################################
#' meanvarsimu
#'
#' @details computes mean and standard deviation of a set of map simulations
#' @param map object generated by genMap
#'
#' @return a vector with 4 elements, the first two give mean and standard deviation of generated raw data, the last two give mean and standard deviation of kriged data
#'
#' @export
#'
#' @examples
#' \donttest{
#' meanvarSimu(mapTest)
#' }
meanvarSimu=function(map)
######################################
{
v=map$rawData
v=v@data[,1]
kv=map$krigData
kv=kv@data[,1]
m=mean(v)
km=mean(kv)
sdd=sd(v)
sdkd=sd(kv)
vec=c(m,km,sdd,sdkd)
names(vec)=c("raw mean","kriged mean","raw sd","kriged sd")
return(vec)
}
###################
#' valZ
#'
#' @details sorts zones according to attribute mean value
#' @param map map object returned by genMap function
#' @param K zoning object (such as returned by calNei function)
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a list with components
#' \describe{
#' \item{val}{list with vector of data values for each zone, zones are sorted by increasing mean values}
#' \item{ord}{order of zones sorted by increasing mean values}
#' }
#'
#' @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)
#' K=criti$zk[[2]][[1]]
#' valZ(mapTest,K)
valZ=function(map,K)
###################
{
tab=map$krigData
pt=K$listZonePoint
mu=K$meanZone
ord=order(mu)
val=list()
k=0
for(ii in ord){
k=k+1
val[[k]]=tab[[1]][pt[[ii]]]
}
return(list(val=val,ord=ord))
}
###################
#' orderZ
#'
#' @details sorts zones according to ord vector
#' @param Z zoning geometry
#' @param ord sorting order
#'
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a zoning geometry (list of SpatialPolygons)
#'
#' @export
#'
#' @examples
#' \donttest{
#' map=genMap(DataObj=NULL,seed=40,disp=FALSE,krig=1,Vnugget=1.2,typeMod="Gau")
#' qProb=c(0.275,0.8)
#' criti=correctionTree(qProb,map,LASTPASS=FALSE)
#' res=searchNODcrit1(qProb,criti)
#' b=res$ind[[1]][1]
#' K=criti$zk[[2]][[b]]
#' Z=K$zonePolygone
#' plotZ(Z)
#' ord=valZ(map,K)$ord
#' Z=orderZ(Z,ord)
#' plotZ(Z)
#' }
###################
orderZ=function(Z,ord)
{
k=0
for (iZ in ord){
k=k+1
Z=setId(Z,iZ,k)
}
return(Z)
}
#############################
#' superLines
#'
#' @details converts boundary (list of x and y pts) into Spatial Lines
#' @param boundary list, contains x and y coordinates of map boundaries
#' @importFrom sp coordinates
#' @importMethodsFrom sp coordinates
#'
#' @return a SpatialLines object
#'
#' @export
#'
#' @examples
#' data(mapTest)
#' superL=superLines(mapTest$boundary)
#' sp::plot(superL)
superLines=function(boundary)
#############################
{
#transform boundary into SpatialLines
boundary = data.frame(boundary)
sp::coordinates(boundary)=~x+y
bl=Line(coordinates(boundary))
bspl=SpatialLines(list(Lines(list(bl),'1')))
bdLines = bspl@lines[[1]]@Lines[[1]]
listBdLines=list(Lines(list(Line(bdLines@coords[1:2,])),'1'))
for (i in 2:(length(bdLines@coords)/2-1))
{
listBdLines[[i]] = Lines(list(Line(bdLines@coords[i:(i+1),])),paste(i))
}
SuperLines = SpatialLines(listBdLines)
return(SuperLines)
}
####################
#' r2
#'
#' @details adjusted R2
#' @param reslm result of a call to lm
#'
#' @return the adjusted r-square of the lm model
#'
#' @export
#' @importFrom stats anova dist lm quantile sd
#'
#' @examples
#' # not run
r2=function(reslm)
####################
{
s2T <- sum(anova(reslm)[[2]]) / sum(anova(reslm)[[1]])
MSE <- anova(reslm)[[3]][2]
adj.R2 <- (s2T - MSE) / s2T
return(adj.R2)
}
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.