如何将netcdf文件中的x、y坐标投影到经纬度上

Question

如何将netcdf文件中的x、y坐标投影到经纬度上

pythonprojectionnetcdfcoordinate-transformationnco

6

我已经从 CCI 网站下载了格陵兰冰盖的速度场 NetCDF 文件。然而，投影方式是这样给出的（如下所示，其中 x 范围为 [-639750,855750]，y 为 [-655750,-3355750]）。如何将这些数据投影到 NetCDF 文件中实际的纬度/经度坐标上？谢谢！如果您感兴趣，可以在这里下载该文件：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'

- Yoni Verhaegen

Python 与此有什么关系？ - WiseDev

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- TomEvans · Accepted Answer

如果您想将整个网格从其本机极射赤面坐标转换为地理（经度和纬度）网格，则可能需要使用类似gdalwarp的工具。不过我认为这不是您要问的问题。

如果我正确理解您的问题，您想从文件中挑选出点，并将它们定位为经度/纬度坐标对。我假设您知道如何从您的netCDF文件中获取一个位置作为XY对，以及该位置的速度值。我还假设您正在使用Python，因为您在此问题上打了这个标签。

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

Pyproj包装了proj4 C库，该库非常广泛地用于坐标系转换。如果您在投影坐标中具有XY对，并且您知道投影坐标系统的定义，则可以像这样使用pyproj的transform函数：

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，并且在此博客文章中有一个示例程序，向您展示如何进行该转换。