pythonnumpyweathermap-projectionsnoaa

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


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: global ocean forecast system

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

Longitude values

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

example


Solution

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

    Regridded data