我在尝试与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相对应，哪个对应于x。可能是其他的，但我不太确定。

Answer 1

总的来说，您的方法对我来说看起来很正确，但是很难看不见您使用的数据，但这里有一些要考虑的要点：

首先，输出文件为空和/或位于错误的位置之间有区别，Georecrencing仅与后者有关。

当互动工作时，您还应确保使用outpation_raster = none正确关闭数据集，这也将为您触发冲洗。

您可以从测试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-island 。

您的方法还假设数据和LAT/LON数组在常规网格上（具有恒定分辨率）。可能并非如此（尤其是如果数据带有坐标的2D网格）。

通常，当给出坐标数组时，它们被定义为对像素的中心有效。与GDAL的Geotransform相比，该晶格定义为像素的（外部）边缘。因此，您可能需要通过减去一半的分辨率来解释这一点。这也会影响您对分辨率的计算，在中心定义的情况下，应使用/（nlons-1）＆amp; /（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)