2

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

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

我的问题是：创建一个 0.5 m x 0.5 m（或 1 x 1 m）的栅格（以 TIF 或 ASCII 格式），其中每个像素的值是 Z 的平均值。如果我在像素中没有点-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

真的感谢所有的帮助

4

1 回答 1

6

一个方法：

  • 定义您的网格spacing（一个浮点数），即同一维度中两个像素/体素中点之间的距离
  • 算出你需要的网格大小，即网格点的数量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)（保存“计数”）
    • 将值添加v到另一个数组中(x_i, y_i)（保存值的总和）
  • 在填充了两个数组后，将值数组之和除以计数数组。0/0.0 将自动分配给NaN，而 c/0.0 将分配给Inf。
于 2012-09-27T12:59:30.757 回答