我有在美国切萨皮克湾不同地点拍摄的物种丰富度的调查数据,我想以图形方式将数据显示为"热图".
我有一个lat/long坐标和丰富度值的数据框,我将其转换为a SpatialPointsDataFrame并使用autoKrige()automap包中的函数生成插值.
首先,任何人都可以评论我是否正确实现了该autoKrige()功能？
其次,我无法绘制数据并覆盖该地区的地图.或者,我可以指定插值网格来反映Bay的边界(如此处所示)吗？关于我如何做到这一点以及我可能从哪里得到这些信息的任何想法？提供网格autoKrige()看起来很容易.
编辑:感谢Paul的超级有用的帖子!这就是我现在拥有的.无法让ggplot接受插值数据和地图投影:
require(rgdal)
require(automap)
#Generate lat/long coordinates and richness data
set.seed(6)
df=data.frame(
lat=sample(seq(36.9,39.3,by=0.01),100,rep=T),
long=sample(seq(-76.5,-76,by=0.01),100,rep=T),
fd=runif(10,0,10))
initial.df=df
#Convert dataframe into SpatialPointsDataFrame
coordinates(df)=~long+lat
#Project latlong coordinates onto an ellipse
proj4string(df)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
#+proj = the type of projection (lat/long)
#+ellps and +datum = the irregularity in the ellipse represented by planet earth
#Transform the projection into Euclidean distances
project_df=spTransform(df, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")) #projInfo(type="proj")
#Perform the interpolation using …我是sqlalchemy的新手.我已经按照教程创建了一个与mysql db关系的现有数据库的自动化
from sqlalchemy import create_engine, MetaData, Column, Table, ForeignKey
from sqlalchemy.ext.automap import automap_base, generate_relationship
from sqlalchemy.orm import relationship, backref
from config import constr, mytables
def _gen_relationship(base, direction, return_fn,
attrname, local_cls, refferred_cls, **kw):
return generate_relationship(base, direction, return_fn, attrname, local_cls, refferred_cls, **kw)
engine = create_engine(constr)
metadata = MetaData()
metadata.reflect(engine, only=mytables)
Base = automap_base(metadata=metadata)
Base.prepare(engine, reflect=True, generate_relationship=_gen_relationship)
Tableclass1 = Base.classes.table1
Tableclass2 = Base.classes.table2
Table2.ID映射到其中一table1列.但是,当我试图使用查询和加入table1和table2,它报告错误,说"无法找到任何外键关系".既然我知道了这两个表的关系,那么在创建类实例之后是否可以声明这种关系？或者有没有办法在查询功能中明确告诉这种关系？谢谢!
我对ggplot比较新,所以如果我的一些问题非常简单或根本无法解决,请原谅我.
我想要做的是生成一个国家的"热图",其中形状的填充是连续的.此外,我有国家的形状.RData.我使用hadley wickham的脚本将我的SpatialPolygon数据转换为数据框.我数据框的long和lat数据现在看起来像这样
head(my_df)
long lat group
6.527187 51.87055 0.1
6.531768 51.87206 0.1
6.541202 51.87656 0.1
6.553331 51.88271 0.1
这个长/拉数据描绘了德国的轮廓.这里省略了数据框的其余部分,因为我认为不需要它.对于某些长/纬度点,我还有第二个数据帧.这看起来像这样
my_fixed_points
long lat value
12.817 48.917 0.04
8.533 52.017 0.034
8.683 50.117 0.02
7.217 49.483 0.0542
我现在想做的是,根据位于该点一定距离内的所有固定点的平均值,为地图的每个点着色.这样我就可以得到整个国家地图的(几乎)连续着色.到目前为止我所拥有的是用ggplot2绘制的国家地图
ggplot(my_df,aes(long,lat)) + geom_polygon(aes(group=group), fill="white") +
geom_path(color="white",aes(group=group)) + coord_equal()
我的第一个想法是生成位于已绘制的地图内的点,然后计算每个生成点的值,my_generated_point如此
value_vector <- subset(my_fixed_points,
spDistsN1(cbind(my_fixed_points$long, my_fixed_points$lat),
c(my_generated_point$long, my_generated_point$lat), longlat=TRUE) < 50,
select = value)
point_value <- mean(value_vector)
我还没有找到一种方法来产生这些点.就像整个问题一样,我甚至不知道是否有可能以这种方式解决.我现在的问题是,是否存在生成这些点的方法和/或是否有另一种方法来解决问题.
解
感谢Paul,我几乎得到了我想要的东西.以下是荷兰样本数据的示例.
library(ggplot2)
library(sp)
library(automap)
library(rgdal)
library(scales)
#get the spatial data …我试图在R中编写一个使用gstat库的代码来创建插值.我已经阅读了gstat手册,基于互联网上的一些例子,我设法编写了这段代码(这只是一部分):
g <- gstat(id="tec", formula=TEC ~ 1, data=data) ##I create an object
v <- variogram(g) # plot the empirical variogram
plot(v)
mod<-vgm(sill=var(data$TEC),model="Sph",range=200,nugget=200) #create the variogram model
v.fit <- fit.variogram(v, model=mod,fit.method=1) #fit the empirical variogram
Theor_variogram=plot(variogram(g),v.fit,main="WLS Model") #plot the theoretical variogram
plot(Theor_variogram)
## Kriging interpolation
p <- predict.gstat(g, model=v.fit, newdata=predGrid)
我的问题是,当我运行最后一个命令(预测)而不是通过普通克里金插值得到结果时,我得到一个反距离加权(IDW).我在gstat手册中读到:"当没有指定变异函数时,反距离加权插值是默认操作.当指定变异函数时,默认预测方法是普通克里金法."
但是,正如您在我的代码中所看到的,我指定了经验和理论变差函数.你知道我为什么一直得到IDW而不是普通克里金吗？它可以与我的坐标类型相关吗？例如,如果我的坐标彼此接近,或者感兴趣的区域太大？任何帮助都非常有用.
在此先感谢迪米特里斯
我使用R功能autokrige的automap包,但我得到了一个错误,我不知道如何解决它.你有什么提示吗？
谢谢!
sp.poidf <- SpatialPointsDataFrame(sp.poi,thresh.df)
proj4string(sp.poidf) <- CRS("+proj=longlat +datum=WGS84")
pro.df=spTransform(sp.poidf, CRS("+proj=merc +zone=32s +datum=WGS84"))
sp.new <- SpatialPoints(new.poi)
proj4string(sp.new) <- CRS("+proj=longlat +datum=WGS84")
pro.new <- spTransform(sp.new, CRS("+proj=merc +zone=32s +datum=WGS84"))
mykri <- autoKrige(mythresh~1,pro.df,newdata=pro.new)
Error in function (classes, fdef, mtable) :
unable to find an inherited method for function "proj4string", for signature "NULL"
我正在尝试使用 automap 包中的 autoKrige() 函数来实现通用克里金法的简单应用。我有一个不规则间隔的测量网格，我想在精细的空间尺度上在它们之间进行插值。示例代码：
library('automap')
# create an irregularly spaced grid
y <-x <-c(-5,-4,-2,-1,-0.5,0,0.5,1,2,4,5)
grid <-expand.grid(x,y)
names(grid) <-c('x', 'y')
# create some measurements, greatest in the centre, with some noise
vals <-apply(grid,1, function(x) {12/(0.1+sqrt(x[1]^2 + x[2]^2))+rnorm(1,2,1.5)})
# get data into sp format
s <-SpatialPointsDataFrame(grid, data.frame(vals))
# make some prediction locations and get them into sp format
pred <-expand.grid(seq(-5,5,by=0.5), seq(-5,5,by=0.5))
pred <-cbind(pred[,1], pred[,2]) # this seems to be needed, not sure why
pred <-SpatialPoints(pred)
# try universal kriging
surf <-autoKrige(vals~x+y, …