如何将 x,y 坐标投影到 netcdf 文件中的 lat/lon

Question

我已经从 CCI 网站下载了格陵兰冰的速度场 sheet 作为 NetCDF 文件。但是，投影如下（见下文，其中 x 的范围在 [-639750,855750] 和 y [-655750,-3355750] 之间）

如何将这些数据投影到 NetCDF 文件中的实际 lat/lon 坐标？已经谢谢了！感兴趣的朋友：可以在这里下载文件：http://products.esa-icesheets-cci.org/products/downloadlist/IV/

Variables:
    crs                                
           Size:       1x1
           Dimensions: 
           Datatype:   int32
           Attributes:
                       grid_mapping_name                     = 'polar_stereographic'
                       standard_parallel                     = 70
                       straight_vertical_longitude_from_pole = -45
                       false_easting                         = 0
                       false_northing                        = 0
                       unit                                  = 'meter'
                       latitude_of_projection_origin         = 90
                       spatial_ref                           = 'PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3413"]]'
    y                                  
           Size:       5401x1
           Dimensions: y
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'Y'
                       long_name     = 'y coordinate of projection'
                       standard_name = 'projection_y_coordinate'
    x                                  
           Size:       2992x1
           Dimensions: x
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'X'
                       long_name     = 'x coordinate of projection'
                       standard_name = 'projection_x_coordinate'

Answer 1

如果您想将整个网格从其原生极地立体坐标转换为地理（经纬度）网格，您可能需要使用像 gdalwarp 这样的工具。不过，我认为这不是你要问的问题。

如果我没看错你的问题，你想从文件中挑选点并将它们定位为 lon/lat 坐标对。我假设您知道如何从 netCDF 文件中获取 XY 对位置以及该位置的速度值。我还假设你在 Python 中这样做，因为你把那个标签放在了这个问题上。

获得 XY 对后，您只需要一个函数（带有一堆参数）将其转换为 lon/lat。您可以在 pyproj 模块中找到该函数。

Pyproj 包装了 proj4 C 库，该库在坐标系转换方面应用非常广泛。如果你在投影坐标中有一个 XY 对并且你知道投影坐标系的定义，你可以像这样使用 pyproj 的转换函数：

import pyproj

# Output coordinates are in WGS 84 longitude and latitude
projOut = pyproj.Proj(init='epsg:4326')

# Input coordinates are in meters on the Polar Stereographic 
# projection given in the netCDF file
projIn = pyproj.Proj('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 
    +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ',
    preserve_units=True)

# here is a coordinate pair near the middle of your data set
x, y = 0.0, -2000000

# transform x,y to lon/lat
lon, lat = pyproj.transform(projIn, projOut, x, y)

# answer: lon = -45.0; lat = 71.6886

...好了。请注意，输出经度为 -45.0，这应该会给您一种温暖的感觉，因为输入 X 坐标为 0，而 -45.0 是数据集投影的中央子午线。如果您想要以弧度而不是度数为单位的答案，请将转换函数中的 radians kwarg 设置为 True.

现在是困难的部分，这实际上是您首先要做的事情 -- 定义用作转换函数参数的 projIn 和 projOut。这些位于转换的输入和输出坐标系中。这些是 Proj 对象，它们包含坐标系变换方程的一堆参数。 proj4 开发人员将它们全部封装在一组整齐的函数中，pyproj 开发人员为它们放置了一个很好的 python 包装器，因此您和我不必跟踪所有细节。在我剩下的所有日子里，我都会感谢他们。

输出坐标系很简单

projOut = pyproj.Proj(init='epsg:4326')

pyproj 库可以从 EPSG 代码构建 Proj 对象。 4326 是 WGS 84 lon/lat 的 EPSG 代码。完成。

设置 projIn 更难，因为您的 netCDF 文件使用 WKT 字符串定义其坐标系，proj4 或 pyproj 无法直接读取该字符串（我很确定）。但是，pyproj.Proj() 将以proj4 参数字符串作为参数。我已经为您提供了此操作所需的操作，因此您可以相信我的

+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

相当于此（直接从您的 netCDF 文件复制）：

PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
  GEOGCS["WGS 84",
    DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Polar_Stereographic"],
  PARAMETER["latitude_of_origin",70],
  PARAMETER["central_meridian",-45],
  PARAMETER["scale_factor",1],
  PARAMETER["false_easting",0],
  PARAMETER["false_northing",0],
  UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
AUTHORITY["EPSG","3413"]]'

如果您希望能够更普遍地执行此操作，您将需要另一个模块来将 WKT 坐标系定义转换为 proj4 参数字符串。 osgeo.osr 就是一个这样的模块，blog post 中有一个示例程序向您展示了如何进行 that 转换。

如何将 x,y 坐标投影到 netcdf 文件中的 lat/lon

How to project x,y coordinates to lat/lon in netcdf file

python

projection

netcdf

coordinate-transformation

nco