将纬度和经度写入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

你可以在Python中使用gdal_grid。我就是这么做的。你只需要像在命令行中那样构建命令，然后把它放在subprocess.call(com, shell=True)里面。记得先导入subprocess模块。

这就是我使用它的方式：

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及以上版本中可用。

Answer 2

虽然这不是你问题中的内容，但看起来你需要将WGS84坐标系的经纬度数据转换成UTM投影。这可以通过GDAL的命令行工具ogr2ogr来实现，使用两个选项-a_srs 4326 -t_srs ????（目标SRID）。你也可以在Python中使用GDAL的OGR模块来完成这个操作。这里有一个使用的例子。

从点数据获取栅格图像有两种独立的方法。第一种是对数据中的值进行插值，这样值就会填满整个区域（有时只填充凸包）。在2D中，有很多方法和工具可以进行插值。使用GDAL时，命令行工具gdal_grid很有用，不过我觉得在Python中可能无法使用。最简单的方式可能是使用scipy.interpolate。一旦你有了一个2D的NumPy数组，使用GDAL/Python来创建栅格文件就很简单了。

第二种将点转换为栅格的方法是将点的位置“烧录”到栅格的像素上。与第一种方法不同，只有点所在的位置有值，而栅格中其他地方的值并没有进行插值。将矢量数据栅格化或烧录到栅格中可以通过GDAL命令行工具gdal_rasterize来完成。你也可以在GDAL/Python中内部完成这个操作，这里有一个例子。