在 MATLAB 中连接地理参考矩阵

Join geo-referenced matrices in MATLAB

我有两个地理参考矩阵(网格化气候数据集),它们具有不同的空间范围和分辨率(即,地球表面一个像素的空间覆盖范围),我想加入。我们称它们为 referenceMatrixtargetMatrix。矩阵可以作为 geotiffs 加载到 MATLAB 中,也可以作为每个像素具有相应 latitudes/longitudes 网格的矩阵加载。

我想要的是用 referenceMatrix 中的相应值填充 targetMatrix 中的 NaN 像素。我可以使用 for 循环逐一查看像素,并根据 referenceMatrix 中的数据填充 targetMatrix 中的 NaNs最近的像素。此处描述了我用来定位 space 中最近像素的方法:
How to find the nearest points to given coordinates with MATLAB?

但是我需要对数千个矩阵执行此操作,而且 for 循环太慢了。因此我想使用逻辑索引,例如

targetMatrix(isnan(targetMatrix)) = referenceMatrix(isnan(targetMatrix))

增加了根据纬度和经度匹配矩阵中像素的能力。你们中的任何人都可以通过基于地理参考比较具有不同范围的矩阵的示例来指出我的方向吗?

下面是输入数据和所需输出的示例

targetMatrix = [1,  NaN, 3; 
                NaN, 5,  6]; 
referenceMatrix = [10, 20, 30, 40; 
                   50, 60, 70, 80]; 
referenceLatitude = [13.3, 13.3, 13.3, 13.3; 
                     14.1, 14.1, 14.1, 14.1]; 
referenceLongitude = [3.2, 4.2, 5.2, 6.2; 
                      3.2, 4.2, 5.2, 6.2];  
targetLatitude = [13.4, 13.4, 13.4; 
                  13.9, 13.9, 13.9]; 
targetLongitude = [3.1, 3.6, 4.1; 
                   3.1, 3.6, 4.1]; 

wantedOutput = [ 1, 10, 3; 
                50,  5, 6];

想要的输出包含来自 targetMatrix 的原始值,其中 NaN 填充了来自 referenceMatrix 的最近(在 space 中)的值,即 1050.

使用isnan. Convert the lat/lon from degrees to radians and then to cartesian coordinates using sph2cart to get actual geodesic distances. Then use knnsearch查找targetMatrix的哪些条目要替换,以找到最接近参考坐标的点的索引相关的目标坐标。使用这些索引从 referenceMatrix 中提取相关条目并用它们替换 NaNs。

nanent = isnan(targetMatrix);
[tarX, tarY] = sph2cart(targetLongitude*pi/180, targetLatitude*pi/180, 1);
[refX, refY] = sph2cart(referenceLongitude*pi/180, referenceLatitude*pi/180, 1);
tmptar = [tarY(nanent) tarX(nanent)];
tmpref = [refY(:) refX(:)];
ind = knnsearch(tmpref, tmptar);
wantedOutput = targetMatrix;
wantedOutput(nanent) = referenceMatrix(ind);

在这种情况下,测试在使用 knnsearch 之前将 lat/lon 转换为笛卡尔坐标 比其他方式加快 knnsearch

knnsearch默认求欧氏距离。您也可以使用 pdist2 and min 的组合来代替 knnsearch(第 6 行)。即

[~, ind] = min(pdist2(tmptar, tmpref), [], 2);

或者您可以在第 6 行中使用 desearchn。即

ind = dsearchn(tmpref, tmptar);

但是 knnsearch 经过测试 在这种情况下比 dsearchn 快。


knnsearchpdist2 需要 ≥ R2010a 以及 Stats 和 ML Toolbox。如果没有,请使用 dsearchn
所有测试均由 OP.

完成