PseudoCode ‐ Geoprocessing

GOES

• Add CRS to xarray datsdet
• Add lat/lon coordinates
• Resample to Uniform Grid
• Interpolate NANs
• Validate

Add CRS to xr.Dataset

We want to add the latitude and longitude coordinates to the SWATH projection.

def add_crs_goes16(ds: xr.Dataset) -> xr.Dataset:
    # get perspective height
    sat_height = ds.goes_imager_projection.attrs["perspective_point_height"]

    # reassign coordinates to correct height
    ds = ds.assign_coords({"x": ds.x.values * sat_height})
    ds = ds.assign_coords({"y": ds.y.values * sat_height})

    # load CRS
    cc = CRS.from_cf(ds.goes_imager_projection.attrs)

    # assign CRS to dataarray
    ds =  ds.rio.write_crs(cc.to_string(), inplace=False)
    
    return ds

Add Lat/Lon Coords

def convert_x_y_to_lat_lon(crs: str, x: list[float], y: list[float]) -> tuple[float, float]:
    transformer = pyproj.Transformer.from_crs(crs, "epsg:4326", always_xy=True)
    lons, lats = transformer.transform(x, y)
    return lons, lats

def convert_lat_lon_to_x_y(crs: str, lon: list[float], lat: list[float]) -> tuple[float, float]:
    transformer = pyproj.Transformer.from_crs("epsg:4326", crs, always_xy=True)
    x, y = transformer.transform(lon, lat
    return x, y

def calc_latlon(ds: xr.Dataset) -> xr.Dataset:
    XX, YY = np.meshgrid(ds.x.data, ds.y.data)
    lons, lats = convert_x_y_to_lat_lon(ds.rio.crs, XX, YY)
    # Check if lons and lons_trans are close in value
    # Set inf to NaN values
    lons[lons == np.inf] = np.nan
    lats[lats == np.inf] = np.nan

    ds = ds.assign_coords({"latitude": (["y", "x"], lats), "longitude": (["y", "x"], lons)})
    ds.latitude.attrs["units"] = "degrees_north"
    ds.longitude.attrs["units"] = "degrees_east"
    return ds

Resample

We want to resample our data from the curvilinear lon-lat grid to a regular lon-lat grid.

import pyinterp.backends.xarray

def create_interp_mesh_2D(ds: xr.DataArray):

    assert len(da.values.shape) == 2

    mesh = pyinterp.RTree()

    # extract coordinates - 2D vector
    X, Y = np.meshgrid(ds.x.values, ds.y.values)

    # ravel them - 1D vectors
    X, Y = X.ravel(), Y.ravel()

    mesh.packing(
        np.vstack((X,Y )).T,
        ds.values.ravel()
    )
    return mesh

def resample_lonlat(da: xr.DataArray, resolution: tuple(float)=(0.25, 0.25)):
    # create interpolation mesh
    mesh = create_interp_mesh(da)

    # create query coordinates
    coords = create_latlon_grid(resolution=resolution)    

    # create meshgrid
    LON, LAT = np.meshgrid(*coords, indexing=“ij”)

    # Inverse Distance Weighting
    idw_eta, neighbors = mesh.inverse_distance_weighting(
        np.vstack((LON.ravel(), LAT.ravel())).T,
        within=True,
        radius=5500,
        k=8,
        num_threads=0,
    )
    idw_eta = idw_eta.reshape(X.shape)

    da_new = xr.Dataset(
        {da.name: (("y", "x"), idw_eta)},
        coords = {"lat": (["y", "x"], lat_coords), 
                  "lon": (["y", "x"], lon_coords)}
    )
    da_new.attrs = da.attrs
    
    return ds_new

Fill NANs

We need to fill what nans we find.

def fillnan_gauss_seidel(da: xr.DataArray, variable: str):
    ds = ds.copy()
    ds["lon"] = ds.lon.assign_attrs(units="degrees_east")
    ds["lat"] = ds.lat.assign_attrs(units="degrees_north")

    da.transpose("lon", "lat", "time")[:, :] = pyinterp.fill.gauss_seidel(
        pyinterp.backends.xarray.Grid3D(da)
    )[1]
    return da

Validation

def transform_360_to_180(coord: np.ndarray) -> np.ndarray:
    return (coord + 180) % 360 - 180

def transform_180_to_90(coord: np.ndarray) -> np.ndarray:
    return (coord + 90) % 180 - 90

def transform_180_to_360(coord: np.ndarray) -> np.ndarray:
    return coord % 360

def transform_90_to_180(coord: np.ndarray) -> np.ndarray:
    return coord % 180

def validate_time(ds: xr.Dataset) -> xr.Dataset:
    ds["time"] = pd.to_datetime(ds.time)
    return ds

def validate_lon_coords(ds: xr.Dataset) -> xr.Dataset:
    new_ds = ds.copy()

    new_ds["lon"] = transform_360_to_180(ds.lon)
    new_ds["lon"] = new_ds.lon.assign_attrs(
        **{
            **dict(
                units="degrees_east",
                standard_name="longitude",
                long_name="Longitude",
            ),
            **ds.lon.attrs,
        }
    )

    return new_ds