在尝试使用GDAL引用图像时，我一直获得空白图像

Question

#!/usr/bin/env python3
import numpy as np
from osgeo import gdal
from osgeo import osr
# Load an array with shape (197, 250, 3)
# Data with dim of 3 contain (value, longitude, latitude)
data = np.load("data.npy")
# Copy the data and coordinates
array = data[:,:,0]
lon = data[:,:,1]
lat = data[:,:,2]
nLons = array.shape[1]
nLats = array.shape[0]
# Calculate the geotransform parameters
maxLon, minLon, maxLat, minLat = [lon.max(), lon.min(), lat.max(), lat.min()]
resLon = (maxLon - minLon) / nLons
resLat = (maxLat - minLat) / nLats
# Get the transform
geotransform = (minLon, resLon, 0, maxLat, 0, -resLat)
# Create the ouptut raster
output_raster = gdal.GetDriverByName('GTiff').Create('myRaster.tif', nLons, nLats, 1, 
gdal.GDT_Int32)
# Set the geotransform
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
# Set to world projection 4326
srs.ImportFromEPSG(4326)
output_raster.SetProjection(srs.ExportToWkt())
output_raster.GetRasterBand(1).WriteArray(array)
output_raster.FlushCache()

上面的代码是为了使用GDAL对栅格进行引用，但是返回空白的tiff文件。我审核数据和变量,,然而,怀疑这个问题可以从geotransform变量。文档要求变量为:

top-left-x, w-e-pixel-resolution, 0,
top-left-y, 0, n-s-pixel-resolution (negative value)

我用了lats和lons不确定哪一个对应x哪一个对应y，可能是别的什么，但我不太确定

Answer 1

总的来说，我认为你的方法是正确的，但如果没有看到你使用的数据，很难判断，但这里有一些要点需要考虑:

首先，输出文件为空和/或位于错误位置之间存在差异，地理引用仅与后者相关。

在交互工作时，您还应该确保使用output_raster = None正确关闭Dataset，这也将为您触发刷新。

您可以从测试GDAL是否读取您想要写入的数据开始。使用如下命令:

ds = gdal.Open('myRaster.tif')
data_from_disk = ds.ReadAsArray()
ds = None
np.testing.assert_array_equal(data_from_disk, array)

如果它们不相同，则可能是数据类型的问题。例如，将接近0的浮点数写成整数，导致它们被裁剪为0，从而显示为"空"。文件。

关于地理参考，您使用的投影具有以度为单位的坐标。如果以弧度为单位，则输出结果接近null岛。

您的方法还假设数据和纬度/经度数组位于规则网格上(具有恒定的分辨率)。但情况可能并非如此(特别是当数据带有二维坐标网格时)。

通常在给定坐标数组时，它们被定义为对像素的中心有效。与GDAL的geotransform相比，它是为像素的(外)边缘定义的。所以你可能需要通过减去一半的分辨率来解释这个问题。这也会影响你对分辨率的计算，在中心分辨率的情况下，应该使用/ (nLons-1)&/ (nLats-1)。或者使用以下命令验证:

# for a regular grid
resLon = lon[0,1] - lon[0,0]
resLat = lat[1,0] - lat[0,0]

当我用一些虚拟数据运行你的代码片段时，它给了我一个正确的输出(忽略上面提到的中心/边缘问题)。

lat, lon = np.mgrid[89:-90:-2, -179:180:2]
array = np.sqrt(lon**2 + lat**2).astype(np.int32)