我已从 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'
最佳答案
如果您想将整个网格从其原生极地立体坐标转换为地理(经纬度)网格,您可能需要使用类似 gdalwarp 的工具。 。不过,我认为这不是你问的问题。
如果我正确地阅读了您的问题,您想从文件中挑选点并将它们定位为经度/纬度坐标对。我假设您知道如何从 netCDF 文件中获取 XY 对位置,以及该位置的速度值。我还假设您是在 Python 中执行此操作,因为您在这个问题上添加了该标签。
一旦你有了一个 XY 对,你只需要一个函数(带有一堆参数)将其转换为经度/纬度。您可以在 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 经度/纬度的 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向您展示如何进行转换。
关于python - 如何将 x,y 坐标投影到 netcdf 文件中的纬度/经度,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/56851980/