在栅格图层上绘制XYZ点

我想在栅格图层上显示一组点，并使它们与栅格的比例相同。我已经看过其他问题，但还没有解决。点值代表30个台站的浓度。栅格图层显示使用浓度生成的插值网格。在下面的代码中，我得到的插值输出覆盖了比我需要映射的更广的区域。我将它剪辑到我用作底图的shp图层上。然后我映射修剪的插值，并添加点。在栅格图层上绘制XYZ点

我已经得到了地图显示的点和栅格，我不清楚如何让点的规模出来相同的栅格。我想我需要以某种方式从插值中提取颜色渐变，并将它们匹配到我的XYZ文件中的值，但我不确定这是否是正确的道路...

我的代码生成目前为止的地图：

# test data for 4 stations 

XYZData<- data.frame(
X = c(577422.4, 579220.9, 580996.0, 584633.1), 
Y = c(4210598, 4211570, 4212198, 4213455),Z = c(7.1,10.2,12,3)) 

library(raster) 
rasterDF <- raster(idw_out) # idw_out is the output interpolation grid 
newproj<-'my projection' # setting projection 
a<-projectRaster(rasterDF, crs=newproj) 

nfullName="./shpdata/this is my base map.shp" 
nlayerName="name of my base map" 
border_rev<-readOGR(dsn=nfullName,layer=nlayerName, verbose=F) # map extent 
x1<-mask(a, border_rev) # clip interpolation to basemap 

h<- colors(); h1<-h[10:100] # attempting to standardize colors used 
p <- plot(x1, col=h1); points(XYZdata, type='p', col=h1)

谢谢！

编辑：这是光栅的总结：

class  : RasterLayer 
dimensions : 317, 425, 134725 (nrow, ncol, ncell) 
resolution : 250, 250 (x, y) 
extent  : 540425, 646675, 4199125, 4278375 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory 
names  : var1.pred 
values  : 0.1547871, 27.05386 (min, max)

这是链接到光栅：dataset

的XYZ数据是这样的：

mydata <- data.frame(Longitude = c(577422.442467061, 579220.875124347, 
+ 580995.982508038, 584633.125414608, 585696.32284728, 589533.313960728, 
+ 587475.548155474, 582963.372460121, 578940.888532248, 582075.111557177 
+), Latitude = c(4210598.27557511, 
+ 4211569.85438234, 4212197.61081488, 4213455.20812117, 4213299.86256043, 
+ 4213030.38550458, 4214284.20259049, 4215457.67119233, 4213298.22263535, 
+ 4217435.01191314 
+), Result = c(7.19, 10.88, 6.59, 5.92, 4.85, 3.83, 4.57, 4.71, 
+ 5.22, 2.07))

任何意见赞赏。

来源

2016-11-05 user1868064

请提供一个可重现的例子 – loki

起初，我创建了一个可重复的例子：

library(rgdal) 
library(raster) 

xyz <- data.frame(x = seq(-160,160, 25), 
        y = seq(-40, 20, 5), 
        z = rnorm(13) 
) 
rasterDF <- raster(ncol = 100, nrow = 50) 
rasterDF[] <- rnorm(ncell(rasterDF))

然后，你必须在XYZ数据转换为SpatialPointsDataFrame（最终spTransform到您'my projection'）：

xyzSP <- SpatialPointsDataFrame(coords = xyz[,1:2], 
           data = xyz, 
           proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

之后，将这些点以下的栅格单元格的数据连接到z值（只是示例的可重现性的解决方法;您的值应该已经相同）

xyzSP$Z <- rasterDF[xyzSP,]

下一步是计算稍后将在图中使用的颜色。因此，你必须正常化与rasterDF的最小值和最大值的Z值：

ncolors <- 100 
xyzSP$color <- terrain.colors(ncolors)[(xyzSP$z - minValue(rasterDF))/ 
             (maxValue(rasterDF)- minValue(rasterDF)) * ncolors]

之后，您可以绘制结果像这样与plot(...., add = TRUE)

plot(rasterDF, col = terrain.colors(ncolors)) 
plot(xyzSP, add = TRUE, bg = xyzSP$color, pch = 21, cex = 1)

加点剧情

来源

2016-11-05 10:35:47 loki

谢谢loki。非常有意义。但是，这不适用于我的实际光栅文件，不知道为什么。当我使用rasterDF（在我的情况下为x1）对z值进行规范化时，出现此错误：terrain.colors（ncolors）中的错误[（xyzSP $ Result - minValue（x1））/（maxValue（x1） - ：仅限于 0可能与负下标混合 – user1868064

你可以提供你的数据（或它的一个样本），以创建一个可重复的例子吗？ – loki

@ loki。我在上面编辑了尽可能多的有关数据的信息。我不知道如何发布栅格的例子，因为它的大。 – user1868064

在栅格图层上绘制XYZ点

回答

相关问题