标签: kriging

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

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

我的问题:

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

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

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

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

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

我有一个名为"seoul032823"的81个观察的每小时PM10数据集.你可以从这里下载.我在这个数据集上进行了普通克里金法,并得到了克里金预测的空间图.我还可以在国家地图上显示观察数据点.但我不能在国家地图上重叠克里格空间预测图.

我想做什么:我想在南韩地图(不是整个韩国)上重叠我的空间预测图.我感兴趣的领域是北纬37.2N到37.7N和经度126.6E到127.2E.这意味着我需要从韩国地图中裁剪这个区域并在此处重叠预测图.我还需要根据浓度值显示原始观测数据点,这些点将遵循空间地图的颜色.例如,我想要这种类型的地图:

我的R代码用于克里金法,并在韩国地图上显示数据点:

library(sp)
library(gstat)
library(automap)
library(rgdal)
library(e1071)
library(dplyr)
library(lattice)

seoul032823 <- read.csv ("seoul032823.csv")

#plotting the pm10 data on Korea Map
library(ggplot2)
library(raster)

seoul032823 <- read.csv ("seoul032823.csv")
skorea<- getData("GADM", country= "KOR", level=1)
plot(skorea)

skorea<- fortify(skorea)
ggplot()+
  geom_map(data= skorea, map= skorea, aes(x=long,y=lat,map_id=id,group=group),
           fill=NA, colour="black") +
  geom_point(data=seoul032823, aes(x=LON, y=LAT), 
             colour= "red", alpha=0.7,na.rm=T) +
  #scale_size(range=c(2,4))+
  labs(title= "PM10 Concentration in Seoul Area at South Korea",
       x="Longitude", y= "Latitude", size="PM10(microgm/m3)")+
  theme(title= element_text(hjust = 0.5,vjust = 1,face= c("bold")))

# Reprojection
coordinates(seoul032823) <- ~LON+LAT …

我正在考虑使用这种方法来插入我所拥有的一些3D点.作为输入,我在一个区域的不同高度处具有大气浓度的气体.我的数据显示为垂直高度每隔几英尺几十英尺的值,但水平分隔数百英尺(所以'列'紧凑的值).

假设在任何给定时间点,值在垂直方向上的变化明显大于在水平方向上的值.

我想用这个假设来执行3D克里金法(作为我可以调整的参数或统计定义的参数 - 或者/或).

我相信scikit学习模块可以做到这一点.如果可以,我的问题是如何创建离散单元输出？也就是说,输出到尺寸为50 x 50 x 1英尺的3D数据网格中.理想情况下,我希望[x_location,y_location,value]的输出与这些(或类似的)距离分开.

不幸的是,我没有太多的时间来玩它,所以我只是希望在深入研究它之前弄清楚这是否可能在Python中.谢谢!

我使用gstat来预测二项式数据,但预测值高于1且低于0.有谁知道我如何处理这个问题？谢谢.

data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x+y
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE

#glm model
glm.lime <- glm(lime~dist+ffreq, meuse, family=binomial(link="logit"))
summary(glm.lime)

#variogram of residuals
var <- variogram(lime~dist+ffreq, data=meuse)
fit.var <- fit.variogram(var, vgm(nugget=0.9, "Sph", range=sqrt(diff(meuse@bbox\[1,\])^2 + diff(meuse@bbox\[2,\])^2)/4, psill=var(glm.lime$residuals)))   
plot(var, fit.var, plot.nu=T)

#universal kriging
kri <- krige(lime~dist+ffreq, meuse, meuse.grid, fit.var)
spplot(kri[1])

lat    long
7.16   124.21
8.6    123.35
8.43   124.28
8.15   125.08

考虑这些坐标,这些坐标对应于测量降雨量数据的气象站.

R中gstat包的介绍使用了meuse数据集.在本教程的某些时候:https://rpubs.com/nabilabd/118172,这些人在这行代码中使用了"meuse.grid":

data("meuse.grid")

我没有这样的文件,我不知道如何创建它,我可以使用这些坐标创建一个吗？或者至少向我指出讨论如何为自定义区域创建自定义网格的材料(即不使用GADM的管理边界).

可能写错了,甚至不知道这个问题是否对R精明的人有意义.不过,我很乐意听到一些方向,或者至少是提示.非常感谢!

R的总瘤和统计数据.

编辑:看到我发布的教程看起来像的样本网格,这是我想做的事情.

编辑2:这种方法是否可行？https://rstudio-pubs-static.s3.amazonaws.com/46259_d328295794034414944deea60552a942.html

我想在ggmap上绘制热图.

  library(ggmap)
  turku<-get_map('turku', zoom=13)
  turkumap<-ggmap(turku, extent="device", legend="topleft")
  turkumap
  turkumap+geom_density2d(mapping=aes(x = lon, y = lat),data = test, )

我得到的错误是:

 Error in (function (x, y, h, n = 25, lims = c(range(x), range(y)))  : 
 bandwidths must be strictly positive

测试变量是:

  test
         lon     lat var1.pred
  1  22.25320 60.4314 -67.04862
  2  22.25332 60.4314 -67.07793
  3  22.25344 60.4314 -67.11007
  4  22.25356 60.4314 -67.14517
  5  22.25368 60.4314 -67.18336
  6  22.25379 60.4314 -67.22478 
  7  22.25391 60.4314 -67.26956
  8  22.25403 60.4314 -67.31783
  9  22.25415 60.4314 -67.36973
  10 22.25427 …

我有以下代码正是我想要的(它是kriging方法的一部分).但问题是它太慢了,我想知道是否有任何选择将for-loop推向numpy？如果我推出numpy.sum,并在那里使用轴参数,它会加速一点点,但显然这不是瓶颈.关于我如何能够将forloop推向numpy以加快速度,还是以其他方式加速它的任何想法？)

# n = 2116
print GRZVV.shape  # (16309, 2116)
print GinvVV.shape  # (2117, 2117) 
VVg = numpy.empty((GRZVV.shape[0]))

for k in xrange(GRZVV.shape[0]):
    GRVV = numpy.empty((n+1, 1))
    GRVV[n, 0] = 1
    GRVV[:n, 0] = GRZVV[k, :]
    EVV = numpy.array(GinvVV * GRVV)  # GinvVV is numpy.matrix
    VVg[k] = numpy.sum(EVV[:n, 0] * VV)

我发布了ndarrays n矩阵的维度来清除一些东西

编辑:VV的形状是2116

什么是一些好的克里金/插值思路/选项,允许高加权点在绘制的R图上轻微加权点？

康涅狄格州有八个县.我找到了质心,想要绘制这8个县的贫困率.其中三个县人口密集(约100万人),其他五个县人口稀少(约10万人).由于这三个人口稠密的县有超过90%的州人口,我希望这三个人口密集的县完全"压倒"地图并影响整个县境内的其他点.

KrigR fields包中的函数有很多参数,也有可以调用的协方差函数,但我不知道从哪里开始？

这里是可重现的代码,可以快速生成硬边界地图,然后是三个不同加权的地图.希望我可以只对这段代码进行更改,但是它可能需要像geoRglm包一样更复杂的东西？三个加权地图中的两个看起来几乎相同,尽管其中一个加权地图是另一个加权地的10倍.

https://raw.githubusercontent.com/davidbrae/swmap/master/20141001%20how%20to%20modify%20the%20Krig%20function%20so%20a%20huge%20weight%20overwhelms%20nearby%20points.R

谢谢!!

编辑:这是我想要的行为的图片示例 -

我每天有61个测量站的降雨量,为期12年(8000平方公里).

目标是创造5公里和25公里分辨率网格日降雨量产品.由于车站的数量很少,并且即使在雨季也不是所有车站都有雨,我选择使用气候变异图.

典型日(irain)的雨量计测量如下,具有由NA表示的少量缺失值.

7.8  4.4 15.4 19.1  5.8    0   42  6.4   21    21     0     0     0  15.6 0     0    10     5   1.2     0  14.4    NA    25  13.2     0   9.2 2   6.6   7.8  13.2  15.4    NA     9     0  15.5     0  18.6     6 0   4.8  10.6     0    9     0  12.4    NA    12     0     3    14 10    10     0    68  21.8    18  14.8   5.4     7     0    NA

作为日常雨量偏斜,变换i相测试立方根变换和日志变换(log1p)为每天单独.然而,对于所有这些日子来说,两种变换都不适合我用shapiro wilk测试测试.因此,我选择正常分数分数变换(NCR),如Grimes&Pardo(2009)降雨量地统计分析所示.并使用了Ashton Shortridge教授的代码

以下代码用于生成季风季节的气候变异函数.请注意,我曾使用超过30%的电台报告下雨的日子.没有找到任何参考.当我达到65%的天数超过阈值时,选择30%.

lag = 3 …

浏览我发现一些工具在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)

任何帮助表示赞赏.