在 R 中切割栅格数据

编程问答 2022-05-13

问题描述

我有一个有趣的案例需要解决，因为我对此类数据没有太多经验，所以我需要一些帮助。

我有大约 1000 个 TIF 文件，每个文件都超过 30MB。这些文件中的每一个都有至少一个我感兴趣的屋顶。我想根据我存储 lon lat 数据的文件只从 tif 中剪切一些特定位置。我设法绘制了 tif 文件，但我正在努力绘制 lon lat 点。作为下一步，我只想切割点附近的区域。可能吗？

为了导入 TIF 文件，我使用了 stars。

import numpy as np
import timeit
import matplotlib.pyplot as plt

def original_code():
    xsize = data_ind.max() + 1
    nodal_values = np.zeros(xsize,dtype=np.float32)
    for nodes in range(xsize):
        nodal_values[nodes] = np.sum(data_values[np.where(data_ind == nodes)[0]])

def much_better():
    _,idx,_ = np.unique(data_ind,return_counts=True,return_inverse=True)
    nodal_values = np.bincount(idx,data_values)

def slightly_better():
    xsize = data_ind.max() + 1
    idx = np.arange(xsize)[:,None] == data_ind
    nodal_values = [np.sum(data_values[idx[i]]) for i in range(xsize)]

sizes = [i*5 for i in range(1,7)]
original_code_times = np.zeros((len(sizes),))
slightly_better_times = np.zeros((len(sizes),))
much_better_times = np.zeros((len(sizes),))
for i,N in enumerate(sizes):
    print(N)
    data_values = np.random.rand(N)
    data_ind = np.random.randint(0,N,N)

    # Divided by 100 repeats to get average
    original_code_times[i] = timeit.timeit(original_code,number=100) / 100
    much_better_times[i] = timeit.timeit(much_better,number=100) / 100
    slightly_better_times[i] = timeit.timeit(slightly_better,number=100) / 100

# Multiply by 1000 to get everything in ms
original_code_times *= 1000
slightly_better_times *= 1000
much_better_times *= 1000

# %%
plt.figure(dpi=120)
plt.title("Small N's")
plt.plot(sizes,original_code_times,label="Original code")
plt.plot(sizes,slightly_better_times,label="Slightly better")
plt.plot(sizes,much_better_times,label="Much better")
plt.ylabel("Time [ms]")
plt.xlabel("N")
plt.xticks(sizes)
plt.legend()
plt.savefig("small_N.png",dpi=120)
plt.show()
plt.close()

以上方法不起作用，知道为什么吗？

绘制数据点后，栅格中是否有仅提取该点区域的功能？

您可以在这里找到文件： https://drive.google.com/drive/folders/1UCgcqCKHQHc5PsbPv95zEAmQ0qwKoYY0?usp=sharing

解决方法

请不要要求我们下载数据。而是使用代码生成一些数据和/或使用 R 附带的文件。例如：

library(terra)
f <- system.file("ex/elev.tif",package="terra")
r <- rast(f)
xy <- spatSample(r,10,xy=TRUE,na.rm=TRUE)
v <- vect(as.matrix(xy[,1:2]),crs=crs(r))

plot(r)
points(v)

在本例中，栅格和点数据具有相同的 CRS。如果您不是这种情况，则应始终变换点（精确且无损），而不是栅格数据。

v <- project(v,crs(r))

我不知道您对“切割数据”的确切含义是什么，但您可以使用 extract 来获取点位置的栅格像元值。

extract(r,v)
#   ID elevation
#1   1       325
#2   2       281
#3   3       324
#4   4       353
#5   5       271
#6   6       306
#7   7       332
#8   8       260
#9   9       270
#10 10       280

你也可以在点周围做一个缓冲区，也许可以使用crop。

geospatialrrr-raster