0

我正在尝试使用 R 执行从伊比利亚半岛收集的数据的插值频率地图。（类似这样的https://gis.stackexchange.com/questions/147660/strange-spatial-interpolation-results-from-ordinary-kriging）

我的问题是，由于 autokrige 函数的属性 new_data 中存在某种错误，该图没有显示插值数据。

https://cran.r-project.org/web/packages/automap/automap.pdf new_data： 包含预测位置的 sp 对象。new_data 可以是点集、网格或多边形。不得包含 NA。如果未提供此对象，则计算默认值。这是通过获取 input_data 的凸包并在该凸包中放置大约 5000 个网格单元来完成的。

我认为问题在于 R 没有很好地读取转换为 poligons 的地图，因为如果我避免使用这个 new_data 属性，我会得到一个 krigging 值的适当图。但我没有获得伊比利亚半岛的良好形状。有人能帮助我吗？我将不胜感激

在这里你可以看到我的数据： http: //pastebin.com/QHjn4qjP

实际代码：现在，由于我将数据坐标转换为 UTM 投影，我没有收到错误消息，但最后一个图没有插值，整个地图显示为一种颜色:(

 setwd("C:/Users/Ruth/Dropbox/R/pruebas")

#Libraries

library(maps)
library(mapdata)
library(automap)
library(sp)
library(maptools)
library(raster)
library(rgdal)

####################MAPA#############

#obtain the grid of the desired country
map<-getData('GADM', country='ESP', level=0)

#convert the grid to projected data
mapa.utm<-spTransform(mapa3,CRSobj =CRS(" +proj=utm +zone=29 +zone=30 +zone=31  +ellps=WGS84"))

###############################Datos#######################

#submit the data
data1<-read.table("FRECUENCIASH.txt",header=T)
head(data1)
attach(data1)
#convert longlat coordinates to UTM
coordinates(data1)<-c("X","Y")
proj4string(data1) = CRS("+proj=longlat +datum=WGS84") 
data1.utm=spTransform(data1, CRS("+proj=utm +zone=29 +zone=30 +zone=31 +ellps=WGS84 "))

######################Kriging interpolation #####################


#Performs an automatic interpolation

kriging_result<-autoKrige(Z~1,data1.utm,mapa.utm,data_variogram = data1.utm) 

#Plot the kriging

result1<-automapPlot(kriging_result$krige_output,"var1.pred",sp.layout = list("sp.points", data1.utm));result1

result2<-plot(kriging_result);result2
4

1 回答 1

1

您得到的错误与您使用未投影数据作为自动映射输入的事实有关。Automap 只能处理投影数据。谷歌搜索map projections应该会给你一些背景信息。要将您的数据投影到合适的投影，您可以使用spTransformsp 包。

它在没有 newdata 的情况下工作的事实是因为在该对象中未设置投影，因此 automap 无法警告您。但是，带有 latlon 输入数据的 automap 的结果并不可靠。

于 2016-04-30T06:55:07.663 回答