python - 将纬度和经度写入 Geotiff 文件

Question

我基本上是想达到这个问题的反面。

我在 WGS84 坐标系中有一组纬度和经度坐标（带有值），我想通过 gdal python 绑定将其写入 geotiff（或只是添加到 gdal 数据集）。

例如，我的起始数据可能是：

lat = np.array([45.345,56.267,23.425])
lon = np.array([134.689,128.774,111.956])
value = np.array([3.0,6.2,2.5])

怎么可能做到这一点？谢谢！

Answer 1

虽然这不是您的问题，但您似乎需要将 WGS84 基准的纬度/经度数据投影到 UTM 投影。这可以使用 GDAL 中的ogr2ogr命令行，使用两个选项-a_srs 4326 -t_srs ????（目标 SRID）。也可以使用 GDAL 的 OGR 模块在 Python 内部完成。这是一个使用示例。

从点数据中获取栅格有两种独立的方法。第一种是对数据中的值进行插值，以便这些值淹没该区域（或有时仅是凸包）。有许多方法和工具可以在 2D 中插值。使用 GDAL，命令行工具gdal_grid可用于此目的，尽管我认为无法从 Python 中使用。可能最简单的方法是使用scipy.interpolate。一旦你有了一个 2D NumPy 数组，用 GDAL/Python创建一个光栅文件就很简单了。

将点转换为栅格的第二种方法是将点位置刻录为栅格上的像素。与第一种方法不同，只有点所在的位置才有值，而值不会在栅格中的其他任何地方插值。可以通过 GDAL 命令行工具gdal_rasterize将矢量栅格化或刻录到栅格中。它也可以在内部使用 GDAL/Python 完成，这里是一个示例。

Answer 2

可以使用 Python 中的 gdal_grid。我正在使用它。您需要做的就是像在命令行中使用它一样构建命令并将其放入 subprocess.call(com, shell=True) 中。您需要先导入子流程模块。

这实际上是我使用它的方式：

pcall= "gdal_grid  --config 'NUM_THREADS=ALL_CPUS GDAL_CACHEMAX=2000'\
       -overwrite -a invdist:power=2.0:smoothing=2.0:radius1=360.0:radius2=360.0\
       -ot UInt16 -of GTiff -outsize %d %d -l %s -zfield 'Z' %s %s "%(npx, npy,\
       lname,ptshapefile,interprasterfile)

subprocess.call(pcall, shell= True)

NUM_THREADS 选项可从 gdal 1.10+