将纬度和经度写入 Geotiff 文件答案

【问题标题】：Write Latitude and Longitude to Geotiff file将纬度和经度写入 Geotiff 文件
【发布时间】：2013-02-06 18:16:40
【问题描述】：

我基本上是想达到与this 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])

如何做到这一点？谢谢！

【问题讨论】：

数据的位置是网格化的还是地图上的点位置？
它们是地图上的随机点位置。

标签： python geolocation tiff gdal

【解决方案1】：

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

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

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

【讨论】：

我发现了一些代码 HERE 将纬度/经度点转换为 UTM 投影点，然后是像素，这可能正在使用 gdal.ApplyGeoTransform 进行一些插值以在给定单元格中分配每个点值？
@shootingstars 你只需要一个有几个像素点的栅格，中间有空白吗？即，栅格化点？数据代表什么？
没错。一旦我“光栅化”了它们，或者将它们放在光栅像素单元中，我就可以使用我在第一条评论中链接到的逆插值方法。与仅使用 gdal_grid 实用程序相比，这是一种糟糕的处理方式吗？
数据是沿船横断面的鲸鱼目击事件，覆盖了北大西洋的大部分地区。
@shootingstars 好的，我很确定您想要插入数据，而不是栅格化它（这些是不同的）。查看更新的答案。

【解决方案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)

从 gdal 1.10+ 开始提供 NUM_THREADS 选项

【讨论】：