python 如何将netcdf文件的投影转换为隆恩和lats的规则网格？

yk9xbfzb 于 2023-04-04 发布在 Python

关注(0)|答案(2)|浏览(199)

我需要做一个插值对象，我输入一个给定的经度和纬度，对象返回最近的海洋表面电流值。我使用的数据集是。你可以通过以下下载最新的预测this link然后点击今天的日期，在底部是一个名为rtofs_glo_uv_YYYYMMDD.tar.gz的文件。如果你解压缩文件，你会得到三个文件，即：

rtofs_glo_2ds_1hrly_uv_20230330_day1.nc
 rtofs_glo_2ds_1hrly_uv_20230330_day2.nc
 rtofs_glo_2ds_1hrly_uv_20230330_day3.nc

然后你可以在python中使用xarray打开它们：

import xarray as xr
from pathlib import Path

download_folder = Path("")

ds = xr.open_mfdataset(download_folder.glob("rtofs*.nc"))

ds
<xarray.Dataset>
Dimensions:     (MT: 27, Y: 3298, X: 4500)
Coordinates:
  * MT          (MT) datetime64[ns] 2023-03-30 ... 2023-04-02
    Longitude   (Y, X) float32 dask.array<chunksize=(3298, 4500), meta=np.ndarray>
    Latitude    (Y, X) float32 dask.array<chunksize=(3298, 4500), meta=np.ndarray>
  * X           (X) int32 1 2 3 4 5 6 7 8 ... 4494 4495 4496 4497 4498 4499 4500
  * Y           (Y) int32 1 2 3 4 5 6 7 8 ... 3292 3293 3294 3295 3296 3297 3298
    Layer       float64 1.0
Data variables:
    u_velocity  (MT, Y, X) float32 dask.array<chunksize=(9, 3298, 4500), meta=np.ndarray>
    v_velocity  (MT, Y, X) float32 dask.array<chunksize=(9, 3298, 4500), meta=np.ndarray>
Attributes:
    CDI:          Climate Data Interface version 1.9.8 (https://mpimet.mpg.de...
    Conventions:  CF-1.0
    history:      Thu Mar 30 09:26:01 2023: cdo merge rtofs_glo_2ds_1hrly_u_v...
    source:       HYCOM archive file
    institution:  National Centers for Environmental Prediction
    title:        HYCOM ATLb2.00
    experiment:   92.8
    CDO:          Climate Data Operators version 1.9.8 (https://mpimet.mpg.de...

这个文件中使用的网格系统与我习惯的非常不同，经度值不是+/-180而是74到1019.12：

ds.Longitude.min().values
array(74.119995, dtype=float32)
ds.Longitude.max().values
array(1019.12, dtype=float32)

ds.Latitude.max().values
array(89.97772, dtype=float32)
ds.Latitude.min().values
array(-78.64, dtype=float32)

有一个different projection being used：

然而，我不确定这些经度值与实际经度之间的关系。
如果我绘制经度值，删除最后10行（因为它们掩盖了细节，使其比其他值大得多），它们看起来像：

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np

ax = plt.subplot()
im = ax.imshow(ds.Longitude.values[:-10, :])
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(im, cax=cax)
plt.show()

我怎样才能改变这个投影，以便我能找到给定经度和纬度的表面电流？
您可以绘制数据集并查看投影：

ds.sel(MT=ds.MT[0]).u_velocity.plot()

python

来源：https://stackoverflow.com/questions/75888712/how-do-i-convert-the-projection-of-a-netcdf-file-to-a-regular-grid-of-lons-and-l

2条答案

按热度按时间

hpxqektj1#

将此数据转换为规则格网的最简单方法是使用CDO
这是一个实用程序函数，用于生成文件，CDO使用该文件定义要插值的网格：

def to_latlon(lon, lat, res):
    if lat[1] < lat[0]:
        raise ValueError("Check lat order")
    if lon[1] < lon[0]:
        raise ValueError("Check lon order")

    with open("gridfile.txt", "w") as f:
        xsize = int((lon[1] - lon[0]) / res[0]) + 1
        ysize = int((lat[1] - lat[0]) / res[1]) + 1
        lon_step = res[0]
        lat_step = res[1]
        f.write("gridtype = lonlat\n")
        f.write("xsize = " + str(xsize) + "\n")
        f.write("ysize = " + str(ysize) + "\n")
        f.write("xfirst = " + str(lon[0]) + "\n")
        f.write("yfirst = " + str(lat[0]) + "\n")
        f.write("xname = " + "lon" + "\n")
        f.write("xlongname = " + "Longitude" + "\n")
        f.write("xunits = " + "degrees_east" + "\n")
        f.write("yname = " + "lat" + "\n")
        f.write("ylongname = " + "Latitude" + "\n")
        f.write("yunits = " + "degrees_north" + "\n")

        f.write("xinc = " + str(lon_step) + "\n")
        f.write("yinc = " + str(lat_step) + "\n")

示例：

to_latlon([-180, 180], [-90, 90], res =[0.083, 0.083])

gridtype = lonlat
xsize = 4338
ysize = 2169
xfirst = -180
yfirst = -90
xname = lon
xlongname = Longitude
xunits = degrees_east
yname = lat
ylongname = Latitude
yunits = degrees_north
xinc = 0.083
yinc = 0.083

然后，您可以重新网格化数据，如下所示：

cdo remapbil,gridfile.txt rtofs_glo_2ds_1hrly_uv_20230330_day3.nc rtofs_glo_2ds_1hrly_uv_20230330_day3_regular.nc

然后你可以用xarray打开这个文件：

path = Path.home() / "Downloads"
files = path.glob("*bil_classic.nc")
ds = xr.open_mfdataset(files)

并绘制：

ds.sel(MT=ds.MT[0], Layer=1).u_velocity.plot()

赞(0）回复(0）举报 2023-04-04

polkgigr2#

您可以使用nctoolkit来完成此操作，它使用CDO作为重新网格化的后端。以下代码应该可以处理此数据：

import nctoolkit as nc
# open the dataset
ds = nc.open_data("foo.nc")
# do the regridding using bilinear default
ds.to_latlon(lon = [-180, 180], lat = [-90, 90], res =[0.083, 0.083])
# view the result
ds.plot()

赞(0）回复(0）举报 2023-04-04

我来回答

python 如何将netcdf文件的投影转换为隆恩和lats的规则网格？

2条答案

相关问题

热门标签

最新问答