按CF-1.0规范生成的NetCDF投影坐标无法被gdalinfo识别

用Python创建NetCDF栅格时，遵循CF-1.0规范（兼容NetCDF GDAL驱动）写入的UTM 19N投影坐标，有时无法被gdalinfo读取识别，进而导致QGIS无法正确投影。该问题与坐标值本身相关：部分坐标值可被gdalinfo正常读取并在QGIS中正确投影，但另一些值则不会被识别为坐标变量。

触发问题的参数 #

width = 7920
height = 7135

transform = rasterio.transform.from_origin(west=541107.0, north=5438536.5, xsize=1.5, ysize=1.5)

x = rasterio.transform.xy(transform, cols=list(range(0, width, 1)), rows=0)[0]
y = rasterio.transform.xy(transform, cols=0, rows=list(range(0, height, 1)))[1]

此时gdalinfo返回结果：

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 7135.0)
Upper Right ( 7920.0,    0.0)
Lower Right ( 7920.0, 7135.0)
Center      ( 3960.0, 3567.5)

可复现问题的示例代码（含正常与异常场景） #

import netCDF4
import numpy as np
import pyproj
import rasterio

fpath = '/D/Data/TEST/test_netcdf_cf.nc'
ds = netCDF4.Dataset(fpath, "w", format="NETCDF4")

# 可被gdalinfo读取的坐标参数
width= 8026
height= 7340

# 无法被gdalinfo读取的坐标参数
width= 7920
height= 7135

# 创建覆盖马尼夸根半岛的仿射变换
# 将像素大小设为1可解决问题
transform = rasterio.transform.from_origin(west=541107.0, north=5438536.5, xsize=1.5, ysize=1.5)

# 使用以下原点坐标也能修复问题...
# 这是怎么回事？GDAL是否要求网格必须规则？
#transform = rasterio.transform.from_origin(west=541107.25, north=5438536.25, xsize=1.5, ysize=1.5)

x = rasterio.transform.xy(transform, cols=list(range(0, width, 1)), rows=0)[0]
y = rasterio.transform.xy(transform, cols=0, rows=list(range(0, height, 1)))[1]

ds.createDimension('y', height)
ds.createDimension('x', width)
y_var = ds.createVariable('y', 'f4', ('y',))
x_var = ds.createVariable('x', 'f4', ('x',))

y_var[:] = y
x_var[:] = x

# 创建CF grid_mapping
grid_mapping = ds.createVariable('grid_mapping', np.int32, ())

grid_mapping.crs_wtk = cf_grid_mapping['crs_wtk']
grid_mapping.semi_major_axis = cf_grid_mapping['semi_major_axis']
grid_mapping.semi_minor_axis = cf_grid_mapping['semi_minor_axis']
grid_mapping.inverse_flattening = cf_grid_mapping['inverse_flattening']
grid_mapping.reference_ellipsoid_name = cf_grid_mapping['reference_ellipsoid_name']
grid_mapping.longitude_of_prime_meridian = cf_grid_mapping['longitude_of_prime_meridian']
grid_mapping.prime_meridian_name = cf_grid_mapping['prime_meridian_name']
grid_mapping.geographic_crs_name = cf_grid_mapping['geographic_crs_name']
grid_mapping.horizontal_datum_name = cf_grid_mapping['horizontal_datum_name']
grid_mapping.projected_crs_name = cf_grid_mapping['projected_crs_name']
grid_mapping.grid_mapping_name = cf_grid_mapping['grid_mapping_name']
grid_mapping.latitude_of_projection_origin = cf_grid_mapping['latitude_of_projection_origin']
grid_mapping.longitude_of_central_meridian = cf_grid_mapping['longitude_of_central_meridian']
grid_mapping.false_easting = cf_grid_mapping['false_easting']
grid_mapping.false_northing = cf_grid_mapping['false_northing']
grid_mapping.scale_factor_at_central_meridian = cf_grid_mapping['scale_factor_at_central_meridian']

y_var.units = 'm'
y_var.standard_name = 'projection_y_coordinate'
y_var.long_name = 'y coordinate of projection'
y_var.axis = 'Y'

x_var.units = 'm'
x_var.standard_name = 'projection_x_coordinate'
y_var.long_name = 'x coordinate of projection'
x_var.axis = 'X'

# 创建波段维度
ds.createDimension('band', 3)
band_var = ds.createVariable('band', 'f4', ('band',), significant_digits=2)

band_var[:] = [350.23, 460.90, 680.7]

band_var.units = 'nm'
band_var.standard_name = 'sensor_band_central_radiation_wavelength'
band_var.axis = 'B'

# 创建数据变量
Rrs_var = ds.createVariable(
    'Rrs',
    'i4',
    ('band', 'y', 'x',),
    fill_value=21474836,
    compression='zlib',
    complevel=1)

Rrs_var.units = 'sr-1'
Rrs_var.standard_name = 'Rrs'
Rrs_var.long_name = 'Remote sensing reflectance'

Rrs_var.grid_mapping = 'grid_mapping'

RrsData = np.random.random(size=(3, height, width))

print(ds.variables['Rrs'].shape)

ds.variables['Rrs'][0:3, 0:height, 0:width] = RrsData

ds.close()

内容的提问来源于stack exchange，提问作者raphidoc