标签: spatial-interpolation

我有在美国切萨皮克湾不同地点拍摄的物种丰富度的调查数据,我想以图形方式将数据显示为"热图".

我有一个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 …

我一直在努力为我的2D矩阵中的"空"像素提供数据.基本上,我理解(但不是很深)插值技术,如反距离加权,克里金,双立方等.我不完全知道起点(在问题陈述或Python案例中).

问题定义: 我有MxN矩阵(规则网格),其中每个像素代表一定的测量值(下图和此图中使用的数据在 这里).我想使用我作为蓝色像素的现有数据插入"问号空间"(白色空间,也包括相同大小但空白的像素)区域的数据.

我的问题:

1)如何插入此数据.任何人都可以给我一个简单的例子(例如3x3矩阵)来清楚地理解这一点吗？

2)有人可以指导我如何在Python环境中执行解决方案的步骤吗？

3)如何使用Python比较精度意义上的插值技术？

4)您是否认为根据数据密度使用不同的插值是个好主意？

我将非常感谢您的回答和建议.

我有一个二维数据数组,有一些缺失值.有三列:

• X
• ÿ
• 强度

我可以在ggplot2中用y绘制x,强度作为色标.

我想平滑颜色之间的过渡,并idw从gstat包中看到了这个功能.idw旨在以二维方式插入NA.它不应该推断,虽然它在技术上确实尊重数据的极限(两个方向都是±20),但它也试图填补绘图边缘的空隙,如下所示:

我想避免在我所拥有的数据限制之外发生任何外推,包括第一张图中所示数据的右下角.

我怎么能实现这个目标？

编辑:这是一个示例数据集.这与上面显示的数据集不完全相同,但它在右下角又包含一个较大的缺失区域.

structure(list(x = c(10L, 15L, -10L, 0L, -5L, -10L, -15L, 0L, 
-15L, 15L, 5L, 10L, -20L, -5L, -15L, -15L, -5L, 5L, 20L, -20L, 
-15L, 20L, -15L, 5L, -5L, -20L, -5L, 15L, 0L, 0L, 15L, 10L, 0L, 
20L, -10L, 5L, 5L, 0L, 20L, 5L, -15L, 5L, -5L, -5L, -15L, -10L, 
-10L, -10L, -5L, -10L, 15L, 20L, 0L, 20L, -15L, 20L, -20L, -15L, 
10L, 15L, 15L, -5L, …

我有许多地理参考水文数据线,每周分辨率:

Station name, Lat, Long, Week 1 average, Week 2 average ... Week 52 average

不幸的是,我也有一些只有月度分辨率的数据:

Station name, Lat, Long, January average, February average ... December average

而不是"重新发明轮子",任何人都可以推荐一个最喜欢的模块,包或技术,以提供每月价值的每周价值的合理插值？线性会很好,但如果我们可以使用坐标来改善基于附近站点的插值,那将会很好.

我用python标记了这篇文章,因为它是我最近使用的语言(尽管不是它的统计函数).如果答案是"使用类似的统计程序r",那就这样吧,但我很好奇python的内容是什么.谢谢!

我对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 …

我的源数据文件是每 15 分钟捕获一次的天气雷达图像。我正在尝试对中间帧进行平滑插值，并尝试从上次测量开始推断运动。

到目前为止，我的方法是从两个连续图像计算密集光流场（cv2.calcOpticalFlowFarneback），然后使用光流和第一张图像插入（cv2.remap）帧。我还使用相同的方法进行了虚拟运动外推。

在示例视频中，您可以看到这种方法的注意事项，在插值期间，插值帧和测量帧之间存在明显的中断，并且在外插部分（从 8 秒开始）完全迷幻。我尝试平滑光学流场和它有一个小的影响，但结果基本保持不变。

我想对接下来要研究的方法有一些指导，并且是否有针对此类问题的任何opencv特定解决方案。

问题一：如何让插值更流畅？问题二：运动时间外推的正确方法是什么？

浏览我发现一些工具在Python中使用克里格网络是pyKriging和高斯过程回归.但是,我无法使它们中的任何一个工作.第一个对我不起作用(甚至不能导入它):

import pyKriging

  File "~/python3.6/site-packages/pyKriging/krige.py", line 142
    except Exception, err: 
                    ^
  SyntaxError: invalid syntax

第二个我不明白如何使用它.我找不到一个简单的工作示例(例如,这个 rroowwllaanndd答案很棒,但遗憾的是数据不再可供下载)

所以我的问题是,如何使用Kriging插入我的数据？我有几个站点数据保存在numpy数组中,如下所示:

2000      1         1         5.0
2000      1         2         3.4
2000      1         3         0.2

列是年 - 月 - 日 - 降水.我有几个这样的数据数组(st1,st2,st3),另一个数组包含每个站的ID和每个站所在的坐标(stid,因此站1位于经度15.6865,北纬62.6420,和等等).

import numpy as np
st1 = np.array([[2000,1,1,5.0],[2000,1,2,3.4],[2000,1,3,0.2]])
st2 = np.array([[2000,1,1,8.2],[2000,1,2,2.5],[2000,1,3,0.0]])
st3 = np.array([[2000,1,1,np.nan],[2000,1,2,4.5],[2000,1,3,1.2]])

stid = np.array([[1,15.6865,62.6420],[2,15.7325,62.1254],[3,16.1035,61.1449]])

我需要的是每天一个数组(或一个3D数组),它包含每天在这样的网格中用Kriging插值的所有站点的数据:

y = np.arange(61,63,0.125)
x = np.arange(14,17,0.125)
X,Y = np.meshgrid(x,y)

任何帮助表示赞赏.

我有一个非均匀定位的图像样本,并希望插入到常规网格,因为(除其他外)大多数图像图形函数需要一个规则的网格.我注意到有一些MatLab函数(例如参见随机像素的图像插值),显然会这样做,但找不到一个R-package.
这是一个简单的例子.

#make up some 2D func
y<-matrix(rep(1:10,10) -.5 + runif(100),nrow=10)
x<-matrix(rep(1:10,10) -.5 + runif(100),nrow=10)
inmat<-sin(x) + cos(y)

因此,值inmat在随机位置上.我想某种outmat<-interpolate(inmat,x,y,gridx,gridy)函数,其中inmat,x和y或者是所有矩阵或全部向量(展开矩阵).

我也看到SciPy有http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html这样做.包中是否有这样的功能,R还是需要从端口SciPy或MatLab代码中移植？

在我使用的show final solution的答案之后编辑了这个问题

我有来自不同来源的非结构化2D数据集,例如: 这些数据集是3 numpy.ndarray(X,Y坐标和Z值).

我的最终目标是在网格上插入这些数据以转换为图像/矩阵.所以,我需要找到插入这些数据的"最佳网格".而且,为此,我需要在该网格的像素之间找到最佳的X和Y步长.

根据点之间的欧氏距离确定步骤:

使用每个点与其最近邻居之间的欧几里德距离的平均值.

• 使用KDTree/ cKDTree来自scipy.spacial用于X,Y数据的构建树.
• 使用query方法k=2获取距离(如果k=1,距离仅为零,因为查询每个点找到自己).

    # Generate KD Tree
    xy = np.c_[x, y]  # X,Y data converted for use with KDTree
    tree = scipy.spacial.cKDTree(xy)  # Create KDtree for X,Y coordinates.

    # Calculate step
    distances, points = tree.query(xy, k=2)  # Query distances for X,Y points
    distances = distances[:, 1:]  # Remove k=1 zero distances
    step = numpy.mean(distances)  # Result

性能调整:

• 使用scipy.spatial.cKDTree而不是scipy.spatial.KDTree因为它真的更快. …

我正在使用opencv重映射功能将图像映射到另一个坐标系.但是,我的初步测试表明插值存在一些问题.在这里,我举一个简单的例子,一个图像的恒定0.1像素移位,到处都是0但位置[50,50].

import cv2
import numpy as np

prvs = np.zeros((100,80), dtype=np.float32)
prvs[50:51, 50:51] = 1.

grid_x, grid_y = np.meshgrid(np.arange(prvs.shape[1]), np.arange(prvs.shape[0]))
grid_y = grid_y.astype(np.float32)
grid_x = grid_x.astype(np.float32) + 0.1

prvs_remapped = cv2.remap(prvs, grid_x, grid_y, interpolation=cv2.INTER_LINEAR)

print(prvs_remapped[50,50])
print(prvs_remapped[50,49])

给

0.90625
0.09375

但是,考虑到线性插值方法,我期望0.9和0.1.我做错了什么还是这个数字问题？有没有更精确的重映射算法？

谢谢.