小编Nem*_*esi的帖子

使用R,我想覆盖一些空间点和多边形,以便为点分配​​我考虑的地理区域的一些属性.

我最常做的是使用命令over的的sp包.我的问题是,我正在处理全球范围内发生的大量地理参考事件,在某些情况下(特别是在沿海地区),经度和纬度组合略微落在国家/地区边界之外.这是一个基于这个非常好的问题的可重复的例子.

## example data
set.seed(1)
library(raster)
library(rgdal)
library(sp)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(0.30*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts1 <- spsample(p2-p, n=3, type="random")
pts2<- spsample(p, n=10, type="random")
pts<-rbind(pts1, pts2)

## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add=TRUE)

# overlay
pts_index<-over(pts, p)

# result
pts_index

#>     ID_1       NAME_1 ID_2           NAME_2 AREA
#>1    NA         <NA> <NA>             <NA>   NA
#>2    NA         <NA> <NA>             <NA>   NA
#>3    NA         <NA> <NA>             <NA>   NA
#>4     1 …

我正在处理.nc文件中组织的气候数据的日常观察.我使用光栅包的stack命令读取它们.每个文件(对应一年)都是具有以下特征的RasterStack元素:

class       : RasterStack 
dimensions  : 360, 720, 259200, 365  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax) 

每一层都是一天价值的栅格.
我想总结各层,以计算月度值.我相信解决方案应该使用calc或stackApply {raster},但我找不到从层x到层y求和的方法或者在求和之前将RasterStack子集化的方法.

我准备了一个只有12层的示例文件(以减小大小).

我不知道如何提出代码,对不起,但它应该是这样的:

library(raster)
setwd("myfolder")
data<-stack(mydata.nc)

datasum<- stackApply(data, ??? ,fun=sum)

谢谢

让我们假设我有一个3d维数组(x,y,z),并希望将我的数据重组为维度矩阵(x*y,z),如:

my_array <- array(1:600, dim=c(10,5,12))
my_matrix<-data.frame()

for (j in 1:5) {
  for (i in 1:10) {
     my_matrix <- rbind (my_matrix, my_array[i,j,1:12])
 }
}

你能建议一种更快更优雅的方式吗？

谢谢

很抱歉打扰你,但是这可能是5个小时我在这个问题上生气,我无法解决这个问题.

我有一个大约37,000个条目的数据集.它们中的每一个都有自己的Lon和Lat坐标值.检查总体值,它们的范围分别如下:纬度(-54.4871,70.66344)和经度(-177.375,178.4419).这绝对合理.

我使用ArcGIS创建了一个包含这37,000个点的shapefile:一切正常.

然后我需要使用R处理这些数据,我用于代码的命令是(maptools包):

cells <- readShapeSpatial('RES',IDvar="id_obj", 
                          proj4string=CRS("+proj=longlat +datum=WGS84"))

但是R给出了一个错误:

validityMethod(as(object,superClass))中的错误:给予不符合数据的地理CRS:2.76663393422e + 145

(我不知道这个数字来自哪里,它不是我数据集的一部分......)

阅读此博客上的其他帖子似乎原因应该是lon或lat的无效数据,但正如我上面提到的,我的数据集不是这种情况.

我尝试创建不同的shapefile,第一个没有投影,使用几个投影(WGS84 Mercator,web mercator ......),但错误始终是相同的...

谢谢你的帮助.

也发布在gis.stackexchange 上

我在 R 的关卡图中绘制 RasterLayer 和 SpatialPolygonsDataFrame 时遇到了一些问题。投影有问题，但我不明白是什么。

在这里您可以找到我的可 重现示例数据和代码：

library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(rasterVis)
library(maptools)

setwd("C:/...path_to/test2")
data<-read.csv("test.csv", header=TRUE)

#creating the raster from a data.frame and giving projection
raster<-rasterFromXYZ(data, crs="+proj=longlat")
raster_proj<-projectRaster(raster, crs="+proj=longlat +datum=WGS84 
+ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

#loading shapes
nuts<-readShapePoly("NUTS_RG_WGS84.shp", proj4string=CRS("+proj=longlat"))
countries<-readShapePoly("countries.shp", proj4string=CRS("+proj=longlat"))

#subsetting the raster
raster_clip<- crop(raster_proj, countries, snap="near")

#plot
p.strip <- list(cex=1.5, lines=1, fontface='bold') 
x.scale <- list(cex=1.5)
y.scale <- list(cex=1.5)
label <- list(labels=list(width=1, cex=1.5), height=0.95)
levelplot(raster_clip, par.settings = RdBuTheme, margin=FALSE,
      at=seq(min(na.omit(values(raster_clip))), max(na.omit(values(raster_clip))), length.out=15),scales=list(x=x.scale, y=y.scale), 
      par.strip.text=p.strip, colorkey=label, …