Python:网格化点 X、Y 和 Z 以提取统计属性

Question

很抱歉这个简单的问题，但我是 Python 的新手，我需要同样的帮助。

我的数据采用点格式:X、Y、Z。其中 X 和 Y 是坐标，z 是值。

我的问题是:创建一个 0.5 m x 0.5 m(或 1 x 1 m)的栅格(TIF 或 ASCII)，其中每个像素的值是 Z 的平均值。如果我没有点在pixel-i 取值为NAN。

我很高兴能得到一些我可以学习和实现的代码的帮助，

提前感谢您的帮助，我真的需要。

我试着研究并写了一段代码:

from osgeo import gdal, osr, ogr
import math, numpy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
import matplotlib.delaunay

tiff = 'test.tif'

gridSize = 0.5
# my area boundary
xmax, xmin = 640000.06, 636999.83
ymax, ymin = 6070000.3, 6066999.86

# number of row and columns
nx = int(math.ceil(abs(xmax - xmin)/gridSize))
ny = int(math.ceil(abs(ymax - ymin)/gridSize))

# Plot the points
plt.scatter(x,y,c=z)
plt.axis([xmin, xmax, ymin, ymax])
plt.colorbar()
plt.show()

# Generate a regular grid.
xi = np.linspace(xmin, xmax, nx)
yi = np.linspace(ymin, ymax, ny)
xi, yi = np.meshgrid(xi, yi)

从这一点来看，我很难理解如何对 x、y、z 点进行索引，以便知道它们下降的位置。我的第一个想法是给网格一个索引并标记点。之后我可以对像素内的点进行平均。空像素(不存在点)是 NAN。

但我不知道这是处理我的数据的正确方式。

之后我写了下面的代码通过GDAL以TIFF格式保存

target_ds = gdal.GetDriverByName('GTiff').Create(tiff, nx,ny, 1, gdal.GDT_Byte) #gdal.GDT_Int32
target_ds.SetGeoTransform((xmin, gridSize, 0,ymax, 0, -gridSize,))

if EPSG == None:
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(EPSG)
    # Make the target raster have the same projection as the source
    target_ds.SetProjection(proj.ExportToWkt())
else:
    # Source has no projection (needs GDAL >= 1.7.0 to work)
    target_ds.SetProjection('LOCAL_CS["arbitrary"]')

target_ds.GetRasterBand(1).WriteArray(numpy.zeros((ny,nx)))
target_ds = None

非常感谢大家的帮助

Answer 1

一个方法:

• 定义您的网格间距(一个 float )，这是同一维度中两个像素/体素中点之间的距离
• 计算出您需要的网格大小，即 x 和 y 维度中的网格点数，N_x 和 N_y。
• 创建两个该大小的 numpy 数组，所有值都为零，例如使用np.zeros([N_x, N_y])
• 遍历你的一组(x，y，v)点和
  • 将每个 (x, y) 对投影到其对应的像素中，通过两个(整数)索引标识:x_i, y_i = tuple([int(c//spacing) for c in (x, y)] )
  • 将 1 添加到位置 (x_i, y_i) 的 one 数组(保留“count”)
  • 将值 v 添加到另一个数组的位置 (x_i, y_i)(保存值的总和)
• 填充两个数组后，将值和数组除以计数数组。 0/0.0 会自动分配给 NaN，而 c/0.0 会自动分配给 Inf。