#' SpatialData to Raster
#' \code{sp2raster}
#' @param sp SpatialPolygon
#' @param mask Raster mask of mesh domain.
#' @param ngrids Number of grid along x direction.
#' @param resolution Resolution, defaul = NULL, resolution = extent / ngrids
#' @param field Index of field
#' @return Raster map
#' @export
sp2raster <- function (sp, mask = get('PIHM.MASK', envir = .pihm),
ngrids=200,
resolution=NULL, field=1) {
if( is.null(mask) ){
ext <- raster::extent (sp)
xlim=ext[1:2]
ylim=ext[3:4]
if ( resolution<=0 || is.null(resolution)){
dx=diff(xlim) / ngrids;
}else{
dx=resolution
}
r <- raster::raster(ext, res=dx)
}else{
r = mask
}
## Rasterize the shapefile
rr <-raster::rasterize(sp, r, field=field)
return(rr)
}
#' Generate the raster mask of Mesh domain
#' \code{PIHM.mask}
#' @param pm \code{PIHM.mesh}
#' @param n Number of grid
#' @param rr Default mask in .pihm environment
#' @param cellsize Resolution, defaul = NULL, resolution = extent / ngrids
#' @param crs Projection parameter
#' @return Raster map
#' @export
PIHM.mask <- function (pm = readmesh(), crs=NULL,
rr = get('PIHM.MASK', envir=.pihm),
n=40000, cellsize=NULL){
# mesh=readmesh(shp=TRUE); ngrids=100; resolution=0
if(is.null(rr)){
spm =sp.mesh2Shape(pm)
sp0=rgeos::gUnaryUnion(spm)
if(is.null(cellsize)){
# grd <- as.data.frame(sp::spsample(spm, "regular", n=n))
# grd <- as.data.frame(sp::spsample(spm, "regular", nsig=2, n=n))
grd <- sp::makegrid(sp0, n = n)
}else{
# grd <- as.data.frame(sp::spsample(spm, "regular", cellsize = cellsize))
grd <- sp::makegrid(sp0, cellsize = cellsize)
}
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE # Create SpatialPixel object
fullgrid(grd) <- TRUE # Create SpatialGrid object
rr=raster::raster(grd); rr[]=1
rr=raster::mask(rr, sp0)
if(!is.null(crs)){
raster::crs(rr) = crs
}
assign('PIHM.MASK', rr, envir=.pihm)
}else{
rr = rr
}
rr
}
#' SpatialData to Raster
#' \code{MeshData2Raster}
#' @param x vector or matrix, length/nrow is number of cells.
#' @param rmask mask of PIHM mesh
#' @param stack Whether export the stack, only when the x is a matrix, i.e. (Ntime x Ncell).
#' @param crs Projejction parameter
#' @param pm PIHM mesh
#' @param plot Whether plot the result.
#' @return Raster map
#' @export
MeshData2Raster <- function(x=getElevation(),
rmask=PIHM.mask(crs=crs),
pm=readmesh(), crs=NULL,
stack=FALSE,
plot =FALSE){
msg = 'MeshData2Raster:'
if(stack){
sl = list()
nx = nrow(x)
for(i in 1:nx){
message(msg, 'RasterStack ', i, '/', nx)
sl[[i]] = PIHMgisR::MeshData2Raster(x[i, ], rmask=rmask, pm=pm,
crs=crs, stack=FALSE, plot=FALSE)
}
ret <- raster::stack(sl)
}else{
if( is.matrix(x) | is.data.frame(x)){
x = as.numeric(x[nrow(x),])
}
if(any(is.na(x)) ){
x[is.na(x)] = 0
}
if (any(is.infinite(x))){
x[is.infinite(x)] = 0
}
xy=getCentroid(pm=pm)[,2:3]
val= data.frame(xy, x)
colnames(val) = c('X', 'Y', 'Z')
sp::coordinates(val) = c('X', 'Y')
grd=methods::as(rmask, 'SpatialGrid')
# if(grepl(method, 'idw')){
# Interpolate the grid cells using a power value of 2 (idp=2.0)
# dat <- gstat::idw(Z ~ 1, val, newdata=grd, idp=2.0)
# r = raster::raster(dat); plot(r)
# }else{
# xlim=diff(range(xy[,1], na.rm = TRUE))
# ylim=diff(range(xy[,2], na.rm = TRUE))
# f.1 <- as.formula(Z ~ 1)
# # Compute the sample variogram; note that the f.1 trend model is one of the
# # parameters passed to variogram(). This tells the function to create the
# # variogram on the de-trended data.
# var.smpl <- gstat::variogram(f.1, val, cloud = FALSE,
# cutoff=xlim/1)
# # Compute the variogram model by passing the nugget, sill and range values
# # to fit.variogram() via the vgm() function.
# dat.fit <- gstat::fit.variogram(var.smpl, fit.ranges = FALSE, fit.sills = FALSE,
# gstat::vgm( model="Gau", nugget=0))
# plot(var.smpl, dat.fit)
# dat <- gstat::krige( f.1, val, grd, dat.fit)
#
# r = raster::raster(dat);plot(r)
# }
tps <- fields::Tps(xy, x)
# use model to predict values at all locations
r <- raster::interpolate(rmask, tps)
# r = raster::raster(dat)
ret <- raster::mask(r,rmask)
}
if(plot){
raster::plot(ret)
}
if(!is.null(crs)){ raster::crs(ret) <- crs }
return(ret)
}
#' Remove the holes in polygons
#' \code{removeholes}
#' @param sp SpatialPolygons or SpatialPolygonDataFrame
#' @return Raster map
#' @export
#' @examples
#' library(sp)
#' p.out = Polygon(cbind(c(4,4,6,7,4),c(5,3,2,5,5)) )
#' p.hole = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4) ), hole = TRUE)
#' sp <- SpatialPolygons(list(Polygons(list(p.out, p.hole), "1")))
#' s = removeholes(sp)
#' par(mfrow=c(1,2))
#' plot(sp)
#' plot(s)
removeholes <- function(sp){
x = sp
nx = length(x)
# vn = rownames(x@data)
rl = list()
for(i in 1:nx){
spg = x@polygons[[i]]@Polygons
npg = length(spg)
ypg = list()
k=1
for(j in 1:npg){
if(spg[[j]]@hole){
}else{
ypg[[k]] = spg[[j]]
k = k+1
}
}
rl[[i]] = sp::Polygons(ypg, ID=i)
}
ysp = sp::SpatialPolygons(rl)
# ret = SpatialPolygonsDataFrame(ysp, data=x@data)
ret = ysp
if( !is.na(raster::crs(x)) & ! is.null(raster::crs(x)) ){
raster::crs(ret) = raster::crs(x)
}
return(ret)
}
#' Generatue fishnet
#' \code{fishnet}
#' @param ext Extension of the fishnet. c(xmin, xmax, ymin, ymax)
#' @param n Number of dx
#' @param dx Interval of x direction
#' @param dy Interval of y direction
#' @param crs Projection
#' @param polygons Whether to export SpatialPolygons
#' @param points Whether to export SpatialPoints
#' @param lines Whether to export SpatialLines
#' @export
#' @examples
#' library(raster)
#' pg=fishnet(ext=c(-80,80, -50,50), dx=5)
#' plot(pg)
fishnet <- function(ext, crs=sp::CRS("+init=epsg:4326"),
n=10,
dx=diff(ext[1:2])/n, dy=dx,
lines=FALSE, polygons=TRUE, points=FALSE){
xmin = ext[1]
xmax = ext[2]
ymin = ext[3]
ymax = ext[4]
dx = min(dx, diff(ext[1:2]))
dy = min(dy, diff(ext[3:4]))
xx=seq(ext[1], ext[2], by=dx)
yy=seq(ext[3], ext[4], by=dy)
nx = length(xx)
ny = length(yy)
if(lines){
vline=cbind(xx, ymin, xx, ymax)
hline = cbind(xmin, yy, xmax, yy)
mat = rbind(vline, hline)
str=paste(paste('(', mat[,1], mat[,2],',', mat[,3], mat[,4], ')'), collapse = ',')
spl=rgeos::readWKT(paste('MULTILINESTRING(', str, ')'))
df = as.data.frame(mat)
colnames(df) = c('x1', 'y1', 'x2','y2')
spdf=sp::SpatialLinesDataFrame(spl, data = df)
ret = spdf
raster::crs(ret) = crs
return(ret)
}
if(polygons){
xy=expand.grid(xx,yy)
xm = matrix(xy[,1], nx,ny)
ym = matrix(xy[,2], nx, ny)
xloc=abind::abind(as.matrix(xm[-nx, -ny]), as.matrix(xm[-nx, -1]), as.matrix(xm[-1, -1]),
as.matrix(xm[-1, -ny]), as.matrix(xm[-nx, -ny]), along=3)
yloc=abind::abind(as.matrix(ym[-nx, -ny]), as.matrix(ym[-nx, -1]), as.matrix(ym[-1, -1]),
as.matrix(ym[-1, -ny]), as.matrix(ym[-nx, -ny]), along=3)
# df=as.data.frame(matrix(0, nrow=(nx-1)*(ny-1), 6))
df=data.frame(as.numeric(apply(xloc, 1:2, min)),
as.numeric(apply(xloc, 1:2, max)),
as.numeric(apply(yloc, 1:2, min)),
as.numeric(apply(yloc, 1:2, max)))
df = data.frame(df, rowMeans(df[,1:2]), rowMeans(df[,1:2+2]) )
colnames(df) = c('xmin','xmax','ymin', 'ymax','xcenter','ycenter')
str=paste('GEOMETRYCOLLECTION(',
paste(paste('POLYGON((',
paste(xm[-nx, -ny], ym[-nx, -ny], ',' ),
paste(xm[-nx, -1], ym[-nx, -1], ','),
paste(xm[-1, -1], ym[-1, -1], ','),
paste(xm[-1, -ny], ym[-1, -ny], ','),
paste(xm[-nx, -ny], ym[-nx, -ny], '' ), '))' )
, collapse =','),
')' )
# str=paste('MULTIPOLYGON(', paste(xt, collapse = ', '), ')')
SRL = rgeos::readWKT(str)
# x1 = x0@polygons[[1]]@Polygons
# SRL =lapply(1:length(x1), function(x, i) {Polygons(list(x[[i]]), ID=i)}, x=x1 )
ret = sp::SpatialPolygonsDataFrame( Sr=SRL, data=df, match.ID = TRUE)
raster::crs(ret) = crs
return(ret)
}
return(NA)
}
#' Add holes into Polygons
#' \code{AddHoleToPolygon}
#' @param poly SpatialPolygons
#' @param hole Hole Polygon
#' @export
AddHoleToPolygon <-function(poly,hole){
# https://stackoverflow.com/questions/29624895/how-to-add-a-hole-to-a-polygon-within-a-spatialpolygonsdataframe
# invert the coordinates for Polygons to flag it as a hole
coordsHole <- hole@polygons[[1]]@Polygons[[1]]@coords
newHole <- sp::Polygon(coordsHole,hole=TRUE)
# punch the hole in the main poly
listPol <- poly@polygons[[1]]@Polygons
listPol[[length(listPol)+1]] <- newHole
punch <-sp::Polygons(listPol,poly@polygons[[1]]@ID)
# make the polygon a SpatialPolygonsDataFrame as the entry
new <- sp::SpatialPolygons(list(punch),proj4string=poly@proj4string)
new <- sp::SpatialPolygonsDataFrame(new,data=as(poly,"data.frame"))
return(new)
}
#' Cut sptialLines with threshold.
#' \code{sp.CutSptialLines}
#' @param sl SpatialLines or SpatialLineDataFrame
#' @param tol Tolerence. If the length of segment is larger than tolerance, cut the segment until the maximum segment is shorter than tolerance.
#' @export
#' @examples
#' library(PIHMgisR)
#' library(sp)
#' x=1:1000/100
#' l1 = Lines(Line(cbind(x, sin(x)) ), ID='a' )
#' sl = SpatialLines( list(l1) )
#' tol1=5;
#' tol2 =2
#' sl1 = sp.CutSptialLines(sl, tol1)
#' sl2 = sp.CutSptialLines(sl, tol2)
#' par(mfrow=c(1,2))
#' plot(sl1, col=1:length(sl1));title(paste0('Tol=', tol1))
#' plot(sl2, col=1:length(sl2));title(paste0('Tol=', tol2))
#'
#' data(sh)
#' riv=sh$riv
#' x = sp.CutSptialLines(riv, tol=5)
#' par(mfrow=c(2,1))
#' plot(riv, col=1:length(riv), lwd=3);
#'
#' plot(riv, col='gray', lwd=3);
#' plot(add=TRUE, x, col=1:length(x))
sp.CutSptialLines <- function(sl, tol){
msg='sp.CutSptialLines::'
ll = rgeos::gLength(sl, byid = TRUE)
if(all(ll < tol) ){
ret = sl
}else{
nsp = length(sl)
xsl = list(); ik=1
for(i in 1:nsp){
sx = sl[i, ]
pxy = extractCoords(sx,unique = TRUE)
np = nrow(pxy)
dacc = cumsum( sp::LineLength(pxy, sum = FALSE))
# dacc =getDist(pxy)
tol = max(c(tol, min(dacc) ) )
len= rgeos::gLength(sx)
if(len > tol){
nsplit = ceiling(len / tol)
}else{
nsplit = 1
}
dd = len / nsplit
v0 = 1 # Vetex 0, Vetex 1
message(msg, i, '/', nsp, '\t', nsplit, '\t', round(dd, 2) )
for(k in 1:nsplit){
if(v0 >=np){
break
}
# message(msg, '\t', k, '/', nsplit)
dk = dd * k
v1 = order(abs(dacc - dk), decreasing = FALSE)[1] + 1
if(v1 + 1>np){
v1 = np
}
message(msg, v0,'\t', v1)
if(v0 == v1){
next;
}
# plot(sl[i, ]);points(pxy); points(pxy[c(v0, v1), ], pch=2, col=2)
xsl[[ik]]= sp::Lines(sp::Line( pxy[c(v0:v1), ]), ID=ik)
ik=ik+1
# points(pxy[v0:v1,], col=k)
v0=v1
}
}
nsl = length(xsl)
tmp = sp::SpatialLines(xsl, proj4string = raster::crs(sl))
ilen = rgeos::gLength(tmp, byid=TRUE)
att=data.frame('INDEX'=1:length(tmp), 'Length'=ilen)
ret = sp::SpatialLinesDataFrame(tmp, data = att)
}
return(ret)
}
#' Extract values on Raster map. The line is a straight line between (0,1).
#' \code{extractRaster}
#' @param r Raster
#' @param xy coordinates of the line, dim=(Npoints, 2); x and y must be in [0, 1]
#' @param ext extension of value xy.
#' @param plot Whether plot result.
#' @importFrom grDevices dev.off graphics.off png rgb topo.colors
#' @importFrom graphics grid hist lines par plot points
#' @importFrom methods as
#' @importFrom stats dist rnorm time
#' @importFrom utils read.table
#' @export
#' @examples
#' library(raster)
# r <- raster(ncol=36, nrow=18)
# r[] <- 1:ncell(r)
# extractRaster(r)
extractRaster<-function(r, xy=NULL, ext = raster::extent(r), plot=T){
if(is.null(xy)){
ndim = dim(r)
x=0:ndim[2] / ndim[2]
y = rep(0.5, length(x))
xy = cbind(x,y)
}
x = ext[1] + xy[,1] * (ext[2]- ext[1] )
y = ext[3] + xy[,2] * (ext[4]- ext[3] )
if(plot){
raster::plot(r);
points(x, y, col=2)
nx=length(x)
points(x,y)
graphics::arrows(x[1], y[1], x[nx], y[nx], lty=3, lwd=1.5, col=2)
# lines(x,y, lwd=1.5, col=2, lty=2)
}
v = raster::extract(r, cbind(x,y))
ret = cbind('x'=x,'y'=y,'z'=v)
return(ret)
}
#' Simplify SpatialData.
#' @param x SpatialData
#' @return Simplified SpatialData
#' @export
SimpleSpatial <-function(x){
# n1=length(x@polygons)
# nj=unlist(lapply(1:n1, function(i){ length(x@polygons[[i]]@Polygons) } ))
# x@polygons[[1]]@Polygons[[1]]@coords
msg='SimpleSpatial'
ni = length(x@polygons)
k=1
sl=list()
for(i in 1:ni){
nj = length(x@polygons[[1]]@Polygons)
for(j in 1:nj){
cd = x@polygons[[i]]@Polygons[[j]]@coords
np=nrow(cd)
message(msg, i,'-',j ,'\t', np)
sl[[k]] = paste('POLYGON((', paste( paste(cd[,1], cd[,2]), collapse = ',' ), '))')
if(k==1){
str = sl[[k]]
}else{
str = paste(str, ',', sl[[k]] )
}
k=k+1
}
}
r=rgeos::readWKT(paste('GEOMETRYCOLLECTION(', str, ')'))
}
#' Simplify SpatialData.
#' \code{PointInDistance}
#' @param pt 2-column coordinates (x,y).
#' @param tol Tolerance
#' @return Index of points that within tol
#' @export
PointInDistance <- function(pt, tol){
msg='PointInDistance'
dm = as.matrix(stats::dist(pt, diag = TRUE))
dm[dm==0]=NA
# View(dm)
dmin=apply(dm, 1, min, na.rm=T)
id=which(dmin < 100)
id1=id2=NULL
tmp1=tmp=dmin[id]
i=2
for(i in 1:length(id)){
if(id[i] %in% id2){next }
id1 = c(id1, i); id1
tmp1[id1] = NA; tmp1
id[which(tmp1 %in% tmp[id1])]
id2=unique(c(id2, id[which(tmp1 %in% tmp[id1])]))
id2
tmp1=tmp
}
cbind('P1'=id[id1], P2=id2)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.