How do I convert the projection of a netcdf file to a regular grid of lons and lats?

Question

I need to make an interpolation object where I enter a given longitude and latitude and the object returns the nearest ocean surface current value. The dataset I am using is . You can download the latest forecast by following this link Then clicking on todays date and at the bottom is a file named rtofs_glo_uv_YYYYMMDD.tar.gz. If you unpack the file, you get three files i.e:

 rtofs_glo_2ds_1hrly_uv_20230330_day1.nc
 rtofs_glo_2ds_1hrly_uv_20230330_day2.nc
 rtofs_glo_2ds_1hrly_uv_20230330_day3.nc

You can then open these in python using 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...

The grid system used in this file is very different to what I am used to, the longitude values are not +/-180 but 74 to 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)

I believe there is a different projection being used:

However I am not sure how these longitude values correlate with the actual longitudes.

If I plot the longitude values, removing the last 10 rows (as they obscure the detail from being much larger than the other values), they look like:

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()

How can I change this projection so that I can find the surface current for a given longitude and latitude?

You can plot the dataset and see the projection as well:

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

Answer 1

The easiest way to convert this data to a regular grid is to use CDO

This is a utility function to generate a file that is used by CDO to define the grid to interpolate for:

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")

Example:

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

You can then regrid the data like:

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

You can then open the file with xarray:

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

and plot:

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

Answer 2

You can do this using nctoolkit, which uses CDO as a backend for regridding. The following code should work with this data:

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()