NOAA-ORR-ERD / gridded

A single API for accessing / working with gridded model results on multiple grid types
https://noaa-orr-erd.github.io/gridded/index.html
The Unlicense
64 stars 14 forks source link

variable.at() for wonky CBOFS SGRID #15

Open acrosby opened 6 years ago

acrosby commented 6 years ago

Hi Guys,

Just started trying this out, and it's great! I am looking at representing CBOFS on a regular grid for a client's tool that expects a regular grid, but of course the CBOFS grid is really weird. Is this something that can be handled elegantly in gridded or am I better off writing a routing to handle this specific case? (Also I couldn't find a bbox funtion so I wrote a quick and dirty one that I am using below.)

Here is an example, where it is obviously picking up on the weird ignored (filled) overlapping parts of the grid:

import numpy as np
import gridded
import matplotlib.pyplot as plt

def get_bbox(grdd):
    if "ugrid" in grdd.grid.info.lower():
        return [nc.grid.node_lon.min(),
                nc.grid.node_lat.min(),
                nc.grid.node_lon.max(),
                nc.grid.node_lat.max()]
    elif "sgrid" in grdd.grid.info.lower():
        return [nc.grid.nodes[:, :, 0].min(),
                nc.grid.nodes[:, :, 1].min(),
                nc.grid.nodes[:, :, 0].max(),
                nc.grid.nodes[:, :, 1].max()]

f = "http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/DBOFS/MODELS/201801/nos.dbofs.fields.f001.20180116.t12z.nc" 
nc = gridded.Dataset(f)
bbox = get_bbox(nc)
ssh = nc.variables["zeta"]
time = ssh.time.get_time_array()
rlon, rlat = np.meshgrid(np.linspace(bbox[0], bbox[2], 500),
                         np.linspace(bbox[1], bbox[3], 500))
rssh = ssh.at(np.vstack((rlon.ravel(),
                         rlat.ravel())).T, 
              time[0]).reshape(rlon.shape)
plt.pcolormesh(rlon, rlat, rssh)

image

acrosby commented 6 years ago

The FMRC appears to be handled a little better, but not perfectly. They both seem to have a grid variable that includes the topology information so I'm not sure why that is.

import numpy as np
import gridded
import matplotlib.pyplot as plt

f = "http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/CBOFS/fmrc/Aggregated_7_day_CBOFS_Fields_Forecast_best.ncd"
nc = gridded.Dataset(f)
bbox = get_bbox(nc)
ssh = nc.variables["zeta"] 
time = ssh.time.get_time_array()
rlon, rlat = np.meshgrid(np.linspace(bbox[0], bbox[2], 500),
                         np.linspace(bbox[1], bbox[3], 500))
rssh = ssh.at(np.vstack((rlon.ravel(),
                         rlat.ravel())).T, 
              time[0]).reshape(rlon.shape)
plt.pcolormesh(rlon, rlat, rssh)

image

ChrisBarker-NOAA commented 6 years ago

@AmyMacFadyen:

Amy, is CBOFS one of those wierdo grids that give us trouble?

@acrosby: a few of the COOPS models have wierd grid topology -- essentially, they have invalid lat-lon for the nodes that are just outside the model domain, which makes it really hard to build the topology properly.

Though I'm not sure that the symptoms you're seeing are due to that...

AmyMacFadyen commented 6 years ago

Yes, CBOFS is a headache. All "non-water" grid points have lat/lon that are just fill values -- they don't make a valid curvilinear grid. If you want to try a different COOPS OFS to see is the problem persists, I think the Tampa Bay one (TBOFS) has a more normal grid.

ChrisBarker-NOAA commented 6 years ago

I suspected as much...

@acrosby: we've thought about filling in those points by extrapolating from the neighbors, but there are edge cases, that can be a challenge. But if you come up with a way to do it, let us know!

acrosby commented 6 years ago

Yes I only see this specific problem with the CBOFS ROMS and I’m sure it’s due to the weird stuff they do with the grid outside of the wet domain. I’ll see if I can come up with a good way using the gridded api that is generally workable if the cbofs grid conditions are inferred or explicitly triggered. On Wed, Jan 17, 2018 at 7:24 PM Amy MacFadyen notifications@github.com wrote:

Yes, CBOFS is a headache. All "non-water" grid points have lat/lon that are just fill values -- they don't make a valid curvilinear grid. If you want to try a different COOPS OFS to see is the problem persists, I think the Tampa Bay one (TBOFS) has a more normal grid.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NOAA-ORR-ERD/gridded/issues/15#issuecomment-358495522, or mute the thread https://github.com/notifications/unsubscribe-auth/AA0zvIIZ-Z7fr5KTAiF7RS0PwfFgyMlGks5tLo8ogaJpZM4RhcIP .

AmyMacFadyen commented 6 years ago

If you do, please let us know! The DBOFS grid has similar issues.

jay-hennen commented 6 years ago

The way I planned to attack this issue was to write a loading routine that basically created a new grid and data arrays from the old one. Iterate through all the points, and put the valid ones into a Nx2 array. Then, use the structured nature of the original grid to create an unstructured grid from only the valid points. Pick only the valid data points in the same way.

rsignell-usgs commented 6 years ago

CBOFS, although a ROMS model, used the Delft3D grid generator, which does not use the same gridding approach as ROMS gridding tools. It allows more flexibility to follow tributaries because it doesn't try to grid the land area, but as Amy said all the lon/lat values on land are NaN.

The right way to interpolate or plot using this grid would be to use the lon_psi and lat_psi points (the corners of the cell), instead of just the lon_rho, lat_rho points (the centers of the cell).

But we don't have too many tools to do this.

(I did try the workaround several years back, trying to extrapolate the lon/lat points adjacent to water cells from the lon/lat values on the water side. This approach fails.