Nothing
R/grid-sp-polygons.R
In icosa: Global Triangular and Penta-Hexagonal Grids Based on Tessellated Icosahedra
Defines functions
oneFace
correct90lat
dateLineBreak
#' Add a \code{\link[sp]{SpatialPolygons}} object to a predefined slot in a \code{\link{trigrid}} or \code{\link{hexagrid}} object
#'
#' @name newsp
#
#' @rdname newsp
#' @param gridObj (\code{\link{trigrid}} or \code{\link{hexagrid}}) An icosahedral grid.
#' @param res (\code{numeric}) The number of points inserted between two vertices, passed to \code{\link{SpPolygons}}.
#'
#' @return A \code{\link{trigrid}} or \code{\link{hexagrid}} object with the new \code{@sp} slot.
#' @examples
#' a<-trigrid(4)
#' a<-newsp(a)
#' plot(a)
#' @exportMethod newsp
setGeneric(
name="newsp",
package="icosa",
def=function(gridObj,res=NULL){
standardGeneric("newsp")
}
)
#' @rdname newsp
setMethod(
"newsp",
signature="trigrid",
definition=function(gridObj, res=NULL){
gridObj@sp<-SpPolygons(gridObj,res=res)
return(gridObj)
}
)
#' Spatial polygons from an icosahedral grid
#'
#' The function will create a \code{\link[sp]{SpatialPolygons}} class 2d representation of the icosahedral grid.
#'
#' @name SpPolygons
#
#' @rdname SpPolygons
#' @param gridObj (\code{\link{trigrid}} or \code{\link{hexagrid}}) An icosahedral grid.
#' @param res (\code{numeric}) The number of points inserted between two vertices, or \code{NULL}, if this is to be set by the package. The default method increases resolution with lower tessellation values, and is higher for higher absolute latitudes.
#' @param ... Arguments passed to class-specific methods.
#'
#' @return A \code{\link[sp]{SpatialPolygons}} class object.
#' @exportMethod SpPolygons
#' @examples
#' a <- trigrid()
#' sp <- SpPolygons(a)
setGeneric(
name="SpPolygons",
def=function(gridObj,...){
standardGeneric("SpPolygons")
}
)
#' @rdname SpPolygons
setMethod(
"SpPolygons",
signature="trigrid",
definition=function(gridObj, res=NULL){
# center back to origin if not there already
if(gridObj@center[1]!=0 | gridObj@center[2]!=0 | gridObj@center[3]!=0){
gridObj<-translate(gridObj,-gridObj@center)
}
zenith <- matrix(c(0,0, gridObj@r), ncol=3, nrow=1)
nadir <- matrix(c(0,0, -gridObj@r), ncol=3, nrow=1)
zenithFace <- locate(gridObj,zenith)
nadirFace <- locate(gridObj,nadir)
# override these if the x86 doesn't find the vertices at the poles
# vertexCoord<-vertices(gridObj, output="polar")
# if(90%in%vertexCoord[,2]){
# zenithFace<-NA
# }
#
# if(-90%in%vertexCoord[,2]){
# nadirFace<-NA
# }
# prepare resolution vector
if(is.null(res)){
# the entire implementations is then
# if(dynamic)
minres <- ceiling(1/prod(gridObj@tessellation)^2*500)
# maxres should be 50 more, whichever it is - doesn't matter for coarse grids,
# and 50 is enough for everything.
maxres <- minres+100
# then the latitudinal correction needs to be added
# the frequency of cells in latitudinal belts
tabBelt <- table(gridObj@belts)
# how many plus vertices are needed in each belt?
plusBelt <- round((maxres-minres)/as.numeric(tabBelt))
# final resolution vector
res <- minres+plusBelt[gridObj@belts]
}else{
if(is.numeric(res)){
res <- rep(res, nrow(gridObj@faces))
}else{
stop("The provided resolution value is not numeric and not 'NULL'. ")
}
}
#extend the faces
v<-gridObj@skeleton$v
f<-gridObj@skeleton$f[as.logical(gridObj@skeleton$aF),1:3]
# res<-30
# there are no pentagons here :)
pent<-0
#extend to make a matrix
temp<- .Call(Cpp_icosa_ExpandBoundariesToCols_, f, v, res, gridObj@center, pent)
#reorder to the outer representation
# temp2<-temp[,gridObj@skeleton$uiF]
put<-gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)]
temp2<-matrix(NA, nrow=nrow(temp), ncol=max(put))
temp2[,put]<-temp
subLog<-!is.na(apply(temp2,2, sum))
temp2<-temp2[,subLog]
# allNames<-paste("F", gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)], sep="")
allNames<-paste("F", put, sep="")
#make a data frame from the matrix
temp2<-data.frame(temp2)
finalList<-lapply(temp2, function(x){
# for(i in 1:length(temp2)){
# x<-temp2[[i]]
#3 lines
faceID<-allNames[x[length(x)]+1]
l<-(length(x)-4)/9
mat<-cbind(
c(x[(0*l+1):(1*l)],x[(3*l+1):(4*l)],x[(6*l+1):(7*l)], x[9*l+1]),
c(x[(1*l+1):(2*l)],x[(4*l+1):(5*l)],x[(7*l+1):(8*l)], x[9*l+2]),
c(x[(2*l+1):(3*l)],x[(5*l+1):(6*l)],x[(8*l+1):(9*l)], x[9*l+3])
)
# after the matrix is properly structured (if res=NULL, then full of weir d things), just omit them
mat<-mat[mat[,1]!= -80000,]
faceMat<-CarToPol(mat,norad=TRUE, origin=gridObj@center)
faceMat<-correct90lat(faceMat)
# do the dateline corrections and some
ofList<-oneFace(faceMat,nadirFace, zenithFace, faceID)
for(j in 1:length(ofList)){
if(ofList[[j]]@hole){
ofList[[j]]@hole<-FALSE
ofList[[j]]@ringDir<-as.integer(1)
ofList[[j]]@coords<-ofList[[j]]@coords[nrow(ofList[[j]]@coords):1,]
}
}
of<-sp::Polygons(ofList, ID=faceID)
# return(of)
# }
})
# plot(NULL, NULL, xlim=c(-180,180), ylim=c(-90,90))
# for(i in 1:length(finalList)){
# # for(i in 430:458){
# tsp<-SpatialPolygons(finalList[i])
# plot(tsp, add=TRUE)
# Sys.sleep(0.1)
# }
endObj<-suppressWarnings(sp::SpatialPolygons(finalList, proj4string=methods::as(gridObj@crs, "CRS")))
return(endObj)
}
)
#' @rdname SpPolygons
setMethod(
"SpPolygons",
signature="hexagrid",
definition=function(gridObj, res=NULL){
# center back to origin if not there already
if(gridObj@center[1]!=0 | gridObj@center[2]!=0 | gridObj@center[3]!=0){
gridObj<-translate(gridObj,-gridObj@center)
}
zenith <- PolToCar(matrix(c(0,90), ncol=2,nrow=1), radius=gridObj@r, origin=gridObj@center)
nadir <- PolToCar(matrix(c(0,-90), ncol=2,nrow=1), radius=gridObj@r, origin=gridObj@center)
zenithFace <- locate(gridObj,zenith)
nadirFace <- locate(gridObj,nadir)
#extend the faces
v<-gridObj@skeleton$v
f<-gridObj@skeleton$vF[as.logical(gridObj@skeleton$aF),]
# pentagons
pent<-sum(is.na(f[,6]))
# based on the outer representation!!!
pentLogInner <- is.na(apply(f, 1, sum))
pentLogOuter <- rep(NA, length(pentLogInner))
pentLogOuter[gridObj@skeleton$aF] <- pentLogInner
# res<-30
# prepare resolution vector
if(is.null(res)){
# the entire implementations is then
# if(dynamic)
minres <- ceiling(1/prod(gridObj@tessellation)^2*500)
# maxres should be 50 more, whichever it is - doesn't matter for coarse grids,
# and 50 is enough for everything.
maxres <- minres+100
# then the latitudinal correction needs to be added
# the frequency of cells in latitudinal belts
# the belts ordered to the outer representation
outerBelt <- gridObj@belts
# innerBelt <- rep(NA, length(outerBelt))
# reorder from outer to inner respresentation
innerBelt <- outerBelt[gridObj@skeleton$aF]
tabBelt <- table(innerBelt)
# how many plus vertices are needed in each belt?
plusBelt <- round((maxres-minres)/as.numeric(tabBelt))
# final resolution vector
res <- minres+plusBelt[innerBelt]
}else{
if(is.numeric(res)){
res <- rep(res, nrow(gridObj@faces))
}else{
stop("The provided resolution value is not numeric and not 'NULL'. ")
}
}
#extend to make a matrix
temp<- .Call(Cpp_icosa_ExpandBoundariesToCols_, f, v, res, gridObj@center, pent)
tempRes <- temp
#reorder to the outer representation
put<-gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)]
temp2<-matrix(NA, nrow=nrow(temp), ncol=max(put))
temp2[,put]<-temp
# subLog<-!is.na(apply(temp2,2, sum))
# need another solution for this - NAs code other things!!!
subLog <- 1:max(put)%in%put
temp2<-temp2[,subLog]
# allNames<-paste("F", gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)], sep="")
allNames<-paste("F", put, sep="")
#make a data frame from the matrix
temp2<-data.frame(temp2)
finalList <- mapply(function(x, pen){
# finalList<-lapply(temp2,function(x){
# finalList<-list()
# for(i in 1:length(temp2)){
# x<-temp2[[i]]
#6 lines
faceID<-allNames[x[length(x)]+1]
l<-(length(x)-4)/(6*3)
# in case of the tessellation=1 object, only pentagons are present!
if(prod(gridObj@tessellation)==1){
l<-(length(x)-4)/(5*3)
}
if((x[length(x)]+1)<(pent+1)){
#if(pen){
#pentagon case
#parts are lines
mat<-cbind(
# x coordinates
c(
x[(0*l+1):(1*l)],
x[(3*l+1):(4*l)],
x[(6*l+1):(7*l)],
x[(9*l+1):(10*l)],
x[(12*l+1):(13*l)]
),
# y coordinates
c(
x[(1*l+1):(2*l)],
x[(4*l+1):(5*l)],
x[(7*l+1):(8*l)],
x[(10*l+1):(11*l)],
x[(13*l+1):(14*l)]
),
# z coordinates
c(
x[(2*l+1):(3*l)],
x[(5*l+1):(6*l)],
x[(8*l+1):(9*l)],
x[(11*l+1):(12*l)],
x[(14*l+1):(15*l)]
)
)
}else{
# hexagon case
#parts are lines
mat<-cbind(
# x coordinates
c(
x[(0*l+1):(1*l)],
x[(3*l+1):(4*l)],
x[(6*l+1):(7*l)],
x[(9*l+1):(10*l)],
x[(12*l+1):(13*l)],
x[(15*l+1):(16*l)]
),
# y coordinates
c(
x[(1*l+1):(2*l)],
x[(4*l+1):(5*l)],
x[(7*l+1):(8*l)],
x[(10*l+1):(11*l)],
x[(13*l+1):(14*l)],
x[(16*l+1):(17*l)]
),
# z coordinates
c(
x[(2*l+1):(3*l)],
x[(5*l+1):(6*l)],
x[(8*l+1):(9*l)],
x[(11*l+1):(12*l)],
x[(14*l+1):(15*l)],
x[(17*l+1):(18*l)])
)
}
# after the matrix is properly structured (if res=NULL, then full of weir d things), just omit them
mat<-mat[mat[,1]!= -80000,]
faceMat<-CarToPol(mat,norad=TRUE, origin=gridObj@center)
faceMat<-correct90lat(faceMat)
# do the dateline corrections and some
ofList<-oneFace(faceMat,nadirFace, zenithFace, faceID)
for(j in 1:length(ofList)){
if(ofList[[j]]@hole){
ofList[[j]]@hole<-FALSE
ofList[[j]]@ringDir<-as.integer(1)
ofList[[j]]@coords<-ofList[[j]]@coords[nrow(ofList[[j]]@coords):1,]
}
}
of<-sp::Polygons(ofList, ID=faceID)
# finalList<-c(finalList, list(of))
return(of)
# }
}, temp2, pentLogOuter)
# })
# plot(NULL, NULL, xlim=c(-180,180), ylim=c(-90,90))
# for(i in 1:length(finalList)){
# # for(i in 430:458){
# tsp<-SpatialPolygons(finalList[i])
# plot(tsp, add=TRUE)
# # Sys.sleep(0.3)
# }
# switch off the warnings
suppressWarnings(nproj <- methods::as(gridObj@crs, "CRS"))
endObj<-sp::SpatialPolygons(finalList, proj4string=nproj)
return(endObj)
}
)
dateLineBreak<-function(oneLine){
#where is the break- if there is any?
endIndex<-which(abs(diff(oneLine[,1]))>350)
if(length(endIndex)>0){
#one line cannot be broken more than once
oneLineSeg1<-Line(oneLine[1:endIndex,,drop=FALSE])
oneLineSeg2<-Line(oneLine[(endIndex+1):nrow(oneLine),,drop=FALSE])
#add the broken lines to the lines of the face
faceLines <- c(list(oneLineSeg1), list(oneLineSeg2))
#if there is no break, do everything normally
}else{
#then to Line object
oneLine <- sp::Line(oneLine)
#add Line to faceLines list
faceLines<-list(oneLine)
}
}
# do something with the 90/-90latitude data - to avoid long0-lat90 every time
correct90lat<-function(faceMat){
if(90%in%abs(faceMat[,2])){
meaning<- which(abs(faceMat[,2])==90)
for(u in 1:length(meaning)){
#average out the one before or after
#take longitudes before to avoid "subscript out of bounds!"
longs<-faceMat[,1]
neighbours<-longs[c(meaning[u]-1,meaning[u]+1)]
faceMat[meaning[u],1]<- mean(neighbours, na.rm=TRUE)
#in case there are no nas
# special case of the icosahedron
if(sum(is.na(neighbours[1:2]))==0){
if(sign(neighbours[1])!=sign(neighbours[2])){
faceMat[meaning[u],1]<-neighbours[1]
}
}
}
}
return(faceMat)
}
oneFace<-function(faceMat, nadirFace, zenithFace,faceID){
#divide the cells that cross the dateline, if there is a difference in long values bigger than 300
boolDiv<-sign(faceMat[,1])==1 & sum(abs(diff(faceMat[,1]))>300)
check<-NULL
#if there is only one deviating sign in the longitude data - assign 180 with the
latSign<-sign(faceMat[,1])
if(sum(latSign==-1)==1){
#and fairly close to 180
if(abs(abs(faceMat[latSign==-1,1])-180)<3){
faceMat[latSign==-1,1]<--180
}
}
if(sum(latSign==1)==1){
if(abs(abs(faceMat[latSign==1,1])-180)<3){
faceMat[latSign==1,1]<-180
}
}
#if the face is polar| that is not a vertex
if(!is.na(nadirFace)){
boolNad<-nadirFace==faceID
}else{
boolNad <- FALSE
}
if(!is.na(zenithFace)){
boolZen<-zenithFace==faceID
}else{
boolZen <- FALSE
}
#in case the selected face is the zenith
if(boolZen){
cutInd<-which(abs(diff(faceMat[,1]))>320)
signChange<-sign(diff(faceMat[,1])[abs(diff(faceMat[,1]))>320])
# filter cases when there is a rounding error
if(length(signChange)>0){
#when positive
if(signChange>0){
insertMat<-matrix(c(-180,90,180,90),ncol=2, byrow=TRUE)
}else{
insertMat<-matrix(c(180,90,-180,90),ncol=2, byrow=TRUE)
}
faceMat<-rbind(faceMat[1:cutInd,],
insertMat,
faceMat[(cutInd+1):nrow(faceMat),]
)
}else{
boolZen <- FALSE
}
}
if(boolNad){
cutInd<-which(abs(diff(faceMat[,1]))>320)
signChange<-sign(diff(faceMat[,1])[abs(diff(faceMat[,1]))>320])
# filter cases when there is a rounding error
if(length(signChange)>0){
#when positive
if(signChange>0){
insertMat<-matrix(c(-180,-90,180,-90),ncol=2, byrow=TRUE)
}else{
insertMat<-matrix(c(180,-90,-180,-90),ncol=2, byrow=TRUE)
}
faceMat<-rbind(faceMat[1:cutInd,],
insertMat,
faceMat[(cutInd+1):nrow(faceMat),]
)
}else{
boolNad <- FALSE
}
}
#if there is a division, and the faces are not polar
if(sum(boolDiv)>0 & sum(boolDiv)<nrow(faceMat) & !boolZen & !boolNad){
cellPart1<-faceMat[boolDiv,,drop=FALSE]
cellPart2<-faceMat[!boolDiv,,drop=FALSE]
#if the divided cells are at 90 deg. lat (vertex on the pole)
if(90%in%abs(faceMat[,2])){
#add the +-180 longitude 90latitude where it is appropriate
#1. first part of the cell
chk<-abs(cellPart1[,2])
if(which(chk==90)==1){
#in case the first row is present
cellPart1<-rbind(c(sign(cellPart1[1,1])*180, cellPart1[1,2]), cellPart1)
}else{
#in case it is the last!
if(which(chk==90)==nrow(cellPart1)){
#put it afterwards
cellPart1<-rbind(cellPart1, c(sign(cellPart1[which(abs(cellPart1[,2])==90),1])*180, cellPart1[which(abs(cellPart1[,2])==90),2]))
}else{
# put the new value after this
cellPart1a<-cellPart1[1:which(chk==90),]
cellPart1b<-rbind(
c(sign(cellPart1[which(abs(cellPart1[,2])==90),1])*180, cellPart1[which(abs(cellPart1[,2])==90),2]),
cellPart1[(which(chk==90)+1):nrow(cellPart1),]
)
cellPart1<-rbind(cellPart1a, cellPart1b)
}
}
#2. second part of the cell
chk <- abs(cellPart2[,2])
# for the icosahedron
if(!90%in%chk) chk<-c(90, chk)
if(which(chk==90)==1){
#in case the first row is present
cellPart2<-rbind(c(sign(cellPart2[1,1])*180, cellPart2[1,2]), cellPart2)
}else{
#put it afterwards
cellPart2<-rbind(cellPart2, c(sign(cellPart2[which(abs(cellPart2[,2])==90),1])*180, cellPart2[which(abs(cellPart2[,2])==90),2]))
}
}
#if either one of these has only one row and it's longitude is 180
if((nrow(cellPart1)==1 & 0.5>abs(abs(cellPart1[1,1])-180)) | (nrow(cellPart2)==1 & 0.5>abs(abs(cellPart2[1,1])-180))){
# aparent divison, but no true divison is required
# turn the sign, and divide not!
faceMat[abs(faceMat[,1])==180,1] <- -faceMat[abs(faceMat[,1])==180,1]
options(warn=-1)
cell<-sp::Polygon(faceMat)
options(warn=0)
oneFace<-list(cell)
# check<-c(check, i)
}else{
#divide regularly
options(warn=-1)
cellPart1<-sp::Polygon(cellPart1)
cellPart2<-sp::Polygon(cellPart2)
options(warn=0)
oneFace<-list(cellPart1, cellPart2)
}
} else{
# if not
options(warn=-1)
cell<-sp::Polygon(faceMat)
options(warn=0)
oneFace<-list(cell)
}
return(oneFace)
}
Try the icosa package in your browser
Any scripts or data that you put into this service are public.
icosa documentation built on March 31, 2023, 8:33 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.
Defines functions oneFace correct90lat dateLineBreak
#' Add a \code{\link[sp]{SpatialPolygons}} object to a predefined slot in a \code{\link{trigrid}} or \code{\link{hexagrid}} object
#'
#' @name newsp
#
#' @rdname newsp
#' @param gridObj (\code{\link{trigrid}} or \code{\link{hexagrid}}) An icosahedral grid.
#' @param res (\code{numeric}) The number of points inserted between two vertices, passed to \code{\link{SpPolygons}}.
#'
#' @return A \code{\link{trigrid}} or \code{\link{hexagrid}} object with the new \code{@sp} slot.
#' @examples
#' a<-trigrid(4)
#' a<-newsp(a)
#' plot(a)
#' @exportMethod newsp
setGeneric(
name="newsp",
package="icosa",
def=function(gridObj,res=NULL){
standardGeneric("newsp")
}
)
#' @rdname newsp
setMethod(
"newsp",
signature="trigrid",
definition=function(gridObj, res=NULL){
gridObj@sp<-SpPolygons(gridObj,res=res)
return(gridObj)
}
)
#' Spatial polygons from an icosahedral grid
#'
#' The function will create a \code{\link[sp]{SpatialPolygons}} class 2d representation of the icosahedral grid.
#'
#' @name SpPolygons
#
#' @rdname SpPolygons
#' @param gridObj (\code{\link{trigrid}} or \code{\link{hexagrid}}) An icosahedral grid.
#' @param res (\code{numeric}) The number of points inserted between two vertices, or \code{NULL}, if this is to be set by the package. The default method increases resolution with lower tessellation values, and is higher for higher absolute latitudes.
#' @param ... Arguments passed to class-specific methods.
#'
#' @return A \code{\link[sp]{SpatialPolygons}} class object.
#' @exportMethod SpPolygons
#' @examples
#' a <- trigrid()
#' sp <- SpPolygons(a)
setGeneric(
name="SpPolygons",
def=function(gridObj,...){
standardGeneric("SpPolygons")
}
)
#' @rdname SpPolygons
setMethod(
"SpPolygons",
signature="trigrid",
definition=function(gridObj, res=NULL){
# center back to origin if not there already
if(gridObj@center[1]!=0 | gridObj@center[2]!=0 | gridObj@center[3]!=0){
gridObj<-translate(gridObj,-gridObj@center)
}
zenith <- matrix(c(0,0, gridObj@r), ncol=3, nrow=1)
nadir <- matrix(c(0,0, -gridObj@r), ncol=3, nrow=1)
zenithFace <- locate(gridObj,zenith)
nadirFace <- locate(gridObj,nadir)
# override these if the x86 doesn't find the vertices at the poles
# vertexCoord<-vertices(gridObj, output="polar")
# if(90%in%vertexCoord[,2]){
# zenithFace<-NA
# }
#
# if(-90%in%vertexCoord[,2]){
# nadirFace<-NA
# }
# prepare resolution vector
if(is.null(res)){
# the entire implementations is then
# if(dynamic)
minres <- ceiling(1/prod(gridObj@tessellation)^2*500)
# maxres should be 50 more, whichever it is - doesn't matter for coarse grids,
# and 50 is enough for everything.
maxres <- minres+100
# then the latitudinal correction needs to be added
# the frequency of cells in latitudinal belts
tabBelt <- table(gridObj@belts)
# how many plus vertices are needed in each belt?
plusBelt <- round((maxres-minres)/as.numeric(tabBelt))
# final resolution vector
res <- minres+plusBelt[gridObj@belts]
}else{
if(is.numeric(res)){
res <- rep(res, nrow(gridObj@faces))
}else{
stop("The provided resolution value is not numeric and not 'NULL'. ")
}
}
#extend the faces
v<-gridObj@skeleton$v
f<-gridObj@skeleton$f[as.logical(gridObj@skeleton$aF),1:3]
# res<-30
# there are no pentagons here :)
pent<-0
#extend to make a matrix
temp<- .Call(Cpp_icosa_ExpandBoundariesToCols_, f, v, res, gridObj@center, pent)
#reorder to the outer representation
# temp2<-temp[,gridObj@skeleton$uiF]
put<-gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)]
temp2<-matrix(NA, nrow=nrow(temp), ncol=max(put))
temp2[,put]<-temp
subLog<-!is.na(apply(temp2,2, sum))
temp2<-temp2[,subLog]
# allNames<-paste("F", gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)], sep="")
allNames<-paste("F", put, sep="")
#make a data frame from the matrix
temp2<-data.frame(temp2)
finalList<-lapply(temp2, function(x){
# for(i in 1:length(temp2)){
# x<-temp2[[i]]
#3 lines
faceID<-allNames[x[length(x)]+1]
l<-(length(x)-4)/9
mat<-cbind(
c(x[(0*l+1):(1*l)],x[(3*l+1):(4*l)],x[(6*l+1):(7*l)], x[9*l+1]),
c(x[(1*l+1):(2*l)],x[(4*l+1):(5*l)],x[(7*l+1):(8*l)], x[9*l+2]),
c(x[(2*l+1):(3*l)],x[(5*l+1):(6*l)],x[(8*l+1):(9*l)], x[9*l+3])
)
# after the matrix is properly structured (if res=NULL, then full of weir d things), just omit them
mat<-mat[mat[,1]!= -80000,]
faceMat<-CarToPol(mat,norad=TRUE, origin=gridObj@center)
faceMat<-correct90lat(faceMat)
# do the dateline corrections and some
ofList<-oneFace(faceMat,nadirFace, zenithFace, faceID)
for(j in 1:length(ofList)){
if(ofList[[j]]@hole){
ofList[[j]]@hole<-FALSE
ofList[[j]]@ringDir<-as.integer(1)
ofList[[j]]@coords<-ofList[[j]]@coords[nrow(ofList[[j]]@coords):1,]
}
}
of<-sp::Polygons(ofList, ID=faceID)
# return(of)
# }
})
# plot(NULL, NULL, xlim=c(-180,180), ylim=c(-90,90))
# for(i in 1:length(finalList)){
# # for(i in 430:458){
# tsp<-SpatialPolygons(finalList[i])
# plot(tsp, add=TRUE)
# Sys.sleep(0.1)
# }
endObj<-suppressWarnings(sp::SpatialPolygons(finalList, proj4string=methods::as(gridObj@crs, "CRS")))
return(endObj)
}
)
#' @rdname SpPolygons
setMethod(
"SpPolygons",
signature="hexagrid",
definition=function(gridObj, res=NULL){
# center back to origin if not there already
if(gridObj@center[1]!=0 | gridObj@center[2]!=0 | gridObj@center[3]!=0){
gridObj<-translate(gridObj,-gridObj@center)
}
zenith <- PolToCar(matrix(c(0,90), ncol=2,nrow=1), radius=gridObj@r, origin=gridObj@center)
nadir <- PolToCar(matrix(c(0,-90), ncol=2,nrow=1), radius=gridObj@r, origin=gridObj@center)
zenithFace <- locate(gridObj,zenith)
nadirFace <- locate(gridObj,nadir)
#extend the faces
v<-gridObj@skeleton$v
f<-gridObj@skeleton$vF[as.logical(gridObj@skeleton$aF),]
# pentagons
pent<-sum(is.na(f[,6]))
# based on the outer representation!!!
pentLogInner <- is.na(apply(f, 1, sum))
pentLogOuter <- rep(NA, length(pentLogInner))
pentLogOuter[gridObj@skeleton$aF] <- pentLogInner
# res<-30
# prepare resolution vector
if(is.null(res)){
# the entire implementations is then
# if(dynamic)
minres <- ceiling(1/prod(gridObj@tessellation)^2*500)
# maxres should be 50 more, whichever it is - doesn't matter for coarse grids,
# and 50 is enough for everything.
maxres <- minres+100
# then the latitudinal correction needs to be added
# the frequency of cells in latitudinal belts
# the belts ordered to the outer representation
outerBelt <- gridObj@belts
# innerBelt <- rep(NA, length(outerBelt))
# reorder from outer to inner respresentation
innerBelt <- outerBelt[gridObj@skeleton$aF]
tabBelt <- table(innerBelt)
# how many plus vertices are needed in each belt?
plusBelt <- round((maxres-minres)/as.numeric(tabBelt))
# final resolution vector
res <- minres+plusBelt[innerBelt]
}else{
if(is.numeric(res)){
res <- rep(res, nrow(gridObj@faces))
}else{
stop("The provided resolution value is not numeric and not 'NULL'. ")
}
}
#extend to make a matrix
temp<- .Call(Cpp_icosa_ExpandBoundariesToCols_, f, v, res, gridObj@center, pent)
tempRes <- temp
#reorder to the outer representation
put<-gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)]
temp2<-matrix(NA, nrow=nrow(temp), ncol=max(put))
temp2[,put]<-temp
# subLog<-!is.na(apply(temp2,2, sum))
# need another solution for this - NAs code other things!!!
subLog <- 1:max(put)%in%put
temp2<-temp2[,subLog]
# allNames<-paste("F", gridObj@skeleton$aF[as.logical(gridObj@skeleton$aF)], sep="")
allNames<-paste("F", put, sep="")
#make a data frame from the matrix
temp2<-data.frame(temp2)
finalList <- mapply(function(x, pen){
# finalList<-lapply(temp2,function(x){
# finalList<-list()
# for(i in 1:length(temp2)){
# x<-temp2[[i]]
#6 lines
faceID<-allNames[x[length(x)]+1]
l<-(length(x)-4)/(6*3)
# in case of the tessellation=1 object, only pentagons are present!
if(prod(gridObj@tessellation)==1){
l<-(length(x)-4)/(5*3)
}
if((x[length(x)]+1)<(pent+1)){
#if(pen){
#pentagon case
#parts are lines
mat<-cbind(
# x coordinates
c(
x[(0*l+1):(1*l)],
x[(3*l+1):(4*l)],
x[(6*l+1):(7*l)],
x[(9*l+1):(10*l)],
x[(12*l+1):(13*l)]
),
# y coordinates
c(
x[(1*l+1):(2*l)],
x[(4*l+1):(5*l)],
x[(7*l+1):(8*l)],
x[(10*l+1):(11*l)],
x[(13*l+1):(14*l)]
),
# z coordinates
c(
x[(2*l+1):(3*l)],
x[(5*l+1):(6*l)],
x[(8*l+1):(9*l)],
x[(11*l+1):(12*l)],
x[(14*l+1):(15*l)]
)
)
}else{
# hexagon case
#parts are lines
mat<-cbind(
# x coordinates
c(
x[(0*l+1):(1*l)],
x[(3*l+1):(4*l)],
x[(6*l+1):(7*l)],
x[(9*l+1):(10*l)],
x[(12*l+1):(13*l)],
x[(15*l+1):(16*l)]
),
# y coordinates
c(
x[(1*l+1):(2*l)],
x[(4*l+1):(5*l)],
x[(7*l+1):(8*l)],
x[(10*l+1):(11*l)],
x[(13*l+1):(14*l)],
x[(16*l+1):(17*l)]
),
# z coordinates
c(
x[(2*l+1):(3*l)],
x[(5*l+1):(6*l)],
x[(8*l+1):(9*l)],
x[(11*l+1):(12*l)],
x[(14*l+1):(15*l)],
x[(17*l+1):(18*l)])
)
}
# after the matrix is properly structured (if res=NULL, then full of weir d things), just omit them
mat<-mat[mat[,1]!= -80000,]
faceMat<-CarToPol(mat,norad=TRUE, origin=gridObj@center)
faceMat<-correct90lat(faceMat)
# do the dateline corrections and some
ofList<-oneFace(faceMat,nadirFace, zenithFace, faceID)
for(j in 1:length(ofList)){
if(ofList[[j]]@hole){
ofList[[j]]@hole<-FALSE
ofList[[j]]@ringDir<-as.integer(1)
ofList[[j]]@coords<-ofList[[j]]@coords[nrow(ofList[[j]]@coords):1,]
}
}
of<-sp::Polygons(ofList, ID=faceID)
# finalList<-c(finalList, list(of))
return(of)
# }
}, temp2, pentLogOuter)
# })
# plot(NULL, NULL, xlim=c(-180,180), ylim=c(-90,90))
# for(i in 1:length(finalList)){
# # for(i in 430:458){
# tsp<-SpatialPolygons(finalList[i])
# plot(tsp, add=TRUE)
# # Sys.sleep(0.3)
# }
# switch off the warnings
suppressWarnings(nproj <- methods::as(gridObj@crs, "CRS"))
endObj<-sp::SpatialPolygons(finalList, proj4string=nproj)
return(endObj)
}
)
dateLineBreak<-function(oneLine){
#where is the break- if there is any?
endIndex<-which(abs(diff(oneLine[,1]))>350)
if(length(endIndex)>0){
#one line cannot be broken more than once
oneLineSeg1<-Line(oneLine[1:endIndex,,drop=FALSE])
oneLineSeg2<-Line(oneLine[(endIndex+1):nrow(oneLine),,drop=FALSE])
#add the broken lines to the lines of the face
faceLines <- c(list(oneLineSeg1), list(oneLineSeg2))
#if there is no break, do everything normally
}else{
#then to Line object
oneLine <- sp::Line(oneLine)
#add Line to faceLines list
faceLines<-list(oneLine)
}
}
# do something with the 90/-90latitude data - to avoid long0-lat90 every time
correct90lat<-function(faceMat){
if(90%in%abs(faceMat[,2])){
meaning<- which(abs(faceMat[,2])==90)
for(u in 1:length(meaning)){
#average out the one before or after
#take longitudes before to avoid "subscript out of bounds!"
longs<-faceMat[,1]
neighbours<-longs[c(meaning[u]-1,meaning[u]+1)]
faceMat[meaning[u],1]<- mean(neighbours, na.rm=TRUE)
#in case there are no nas
# special case of the icosahedron
if(sum(is.na(neighbours[1:2]))==0){
if(sign(neighbours[1])!=sign(neighbours[2])){
faceMat[meaning[u],1]<-neighbours[1]
}
}
}
}
return(faceMat)
}
oneFace<-function(faceMat, nadirFace, zenithFace,faceID){
#divide the cells that cross the dateline, if there is a difference in long values bigger than 300
boolDiv<-sign(faceMat[,1])==1 & sum(abs(diff(faceMat[,1]))>300)
check<-NULL
#if there is only one deviating sign in the longitude data - assign 180 with the
latSign<-sign(faceMat[,1])
if(sum(latSign==-1)==1){
#and fairly close to 180
if(abs(abs(faceMat[latSign==-1,1])-180)<3){
faceMat[latSign==-1,1]<--180
}
}
if(sum(latSign==1)==1){
if(abs(abs(faceMat[latSign==1,1])-180)<3){
faceMat[latSign==1,1]<-180
}
}
#if the face is polar| that is not a vertex
if(!is.na(nadirFace)){
boolNad<-nadirFace==faceID
}else{
boolNad <- FALSE
}
if(!is.na(zenithFace)){
boolZen<-zenithFace==faceID
}else{
boolZen <- FALSE
}
#in case the selected face is the zenith
if(boolZen){
cutInd<-which(abs(diff(faceMat[,1]))>320)
signChange<-sign(diff(faceMat[,1])[abs(diff(faceMat[,1]))>320])
# filter cases when there is a rounding error
if(length(signChange)>0){
#when positive
if(signChange>0){
insertMat<-matrix(c(-180,90,180,90),ncol=2, byrow=TRUE)
}else{
insertMat<-matrix(c(180,90,-180,90),ncol=2, byrow=TRUE)
}
faceMat<-rbind(faceMat[1:cutInd,],
insertMat,
faceMat[(cutInd+1):nrow(faceMat),]
)
}else{
boolZen <- FALSE
}
}
if(boolNad){
cutInd<-which(abs(diff(faceMat[,1]))>320)
signChange<-sign(diff(faceMat[,1])[abs(diff(faceMat[,1]))>320])
# filter cases when there is a rounding error
if(length(signChange)>0){
#when positive
if(signChange>0){
insertMat<-matrix(c(-180,-90,180,-90),ncol=2, byrow=TRUE)
}else{
insertMat<-matrix(c(180,-90,-180,-90),ncol=2, byrow=TRUE)
}
faceMat<-rbind(faceMat[1:cutInd,],
insertMat,
faceMat[(cutInd+1):nrow(faceMat),]
)
}else{
boolNad <- FALSE
}
}
#if there is a division, and the faces are not polar
if(sum(boolDiv)>0 & sum(boolDiv)<nrow(faceMat) & !boolZen & !boolNad){
cellPart1<-faceMat[boolDiv,,drop=FALSE]
cellPart2<-faceMat[!boolDiv,,drop=FALSE]
#if the divided cells are at 90 deg. lat (vertex on the pole)
if(90%in%abs(faceMat[,2])){
#add the +-180 longitude 90latitude where it is appropriate
#1. first part of the cell
chk<-abs(cellPart1[,2])
if(which(chk==90)==1){
#in case the first row is present
cellPart1<-rbind(c(sign(cellPart1[1,1])*180, cellPart1[1,2]), cellPart1)
}else{
#in case it is the last!
if(which(chk==90)==nrow(cellPart1)){
#put it afterwards
cellPart1<-rbind(cellPart1, c(sign(cellPart1[which(abs(cellPart1[,2])==90),1])*180, cellPart1[which(abs(cellPart1[,2])==90),2]))
}else{
# put the new value after this
cellPart1a<-cellPart1[1:which(chk==90),]
cellPart1b<-rbind(
c(sign(cellPart1[which(abs(cellPart1[,2])==90),1])*180, cellPart1[which(abs(cellPart1[,2])==90),2]),
cellPart1[(which(chk==90)+1):nrow(cellPart1),]
)
cellPart1<-rbind(cellPart1a, cellPart1b)
}
}
#2. second part of the cell
chk <- abs(cellPart2[,2])
# for the icosahedron
if(!90%in%chk) chk<-c(90, chk)
if(which(chk==90)==1){
#in case the first row is present
cellPart2<-rbind(c(sign(cellPart2[1,1])*180, cellPart2[1,2]), cellPart2)
}else{
#put it afterwards
cellPart2<-rbind(cellPart2, c(sign(cellPart2[which(abs(cellPart2[,2])==90),1])*180, cellPart2[which(abs(cellPart2[,2])==90),2]))
}
}
#if either one of these has only one row and it's longitude is 180
if((nrow(cellPart1)==1 & 0.5>abs(abs(cellPart1[1,1])-180)) | (nrow(cellPart2)==1 & 0.5>abs(abs(cellPart2[1,1])-180))){
# aparent divison, but no true divison is required
# turn the sign, and divide not!
faceMat[abs(faceMat[,1])==180,1] <- -faceMat[abs(faceMat[,1])==180,1]
options(warn=-1)
cell<-sp::Polygon(faceMat)
options(warn=0)
oneFace<-list(cell)
# check<-c(check, i)
}else{
#divide regularly
options(warn=-1)
cellPart1<-sp::Polygon(cellPart1)
cellPart2<-sp::Polygon(cellPart2)
options(warn=0)
oneFace<-list(cellPart1, cellPart2)
}
} else{
# if not
options(warn=-1)
cell<-sp::Polygon(faceMat)
options(warn=0)
oneFace<-list(cell)
}
return(oneFace)
}
Try the icosa 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.