Accessing grid coordinates

[kwilcox]

@ayan-usgs @rsignell-usgs In sci-wms, we need to create a tree index for each grid in a sgrid dataset so we can quickly identify the nearest grid point for a given lat/lon pair. When a request comes in for a variable, we look to see which grid it is on (rho, psi, u, or v for ROMS, but it should be arbitrary) and use the appropriate tree index to query for the nearest grid point.

Does access to this information exist in pysgrid? I couldn't find a way to pull it out.

What I'm looking to do:

In [1]: sg.grids
Out[1]: ['rho', 'psi', 'u', 'v']

In [2]: sg.grids['rho'].dimensions
Out[2]: ['eta_rho', 'xi_rho']

In [3]: sg.grids['rho'].variables
Out[3]: ['lat_rho', 'lon_rho']

In [4]: sg.grids['rho'].coordinates
Out[4]: [
            [
                ... lat_rho ...
            ],
            [
                ... lon_rho ...
            ]
        ]

In [5]: sg.grids['rho'].coordinates.x
Out[5]: [
            ... lon_rho ...
        ]

In [6]: sg.grids['rho'].coordinates.y
Out[6]: [
            ... lat_rho ...
        ]

In [7] sg.temperature.grid    # Returns same object as sg.grids['rho']
Out[7] <pysgrid.grid.Grid at 0x7f190ef7d950>

In [8] sg.temperature.grid.dimensions
Out[8]: ['eta_rho', 'xi_rho']

In [9]:  sg.temperature.grid.coordinates
Out[9]: [
            [
                ... lat_rho ...
            ],
            [
                ... lon_rho ...
            ]
        ]

What is .centers actually returning? Which grids? All that are available? How can I choose which to get back?

I think this is close to what I'm looking for

In [4]: sg.centers.T
Out[4]: 
array([[[ -93.51105439,  -93.53564794,  -93.5602415 , ..., -101.70070766,
         -101.72530121, -101.74989477],
        [ -93.46607342,  -93.49066489,  -93.51525636, ..., -101.65503347,
         -101.67962494, -101.70421642],
        [ -93.42109245,  -93.44568184,  -93.47027123, ..., -101.60935928,
         -101.63394867, -101.65853806],
        ..., 
        [ -53.34304631,  -53.36578088,  -53.38851544, ...,  -60.91365773,
          -60.9363923 ,  -60.95912687],
        [ -53.29806533,  -53.32079782,  -53.34353031, ...,  -60.86798354,
          -60.89071603,  -60.91344852],
        [ -53.25308436,  -53.27581477,  -53.29854517, ...,  -60.82230936,
          -60.84503976,  -60.86777017]],

       [[  11.88846735,   11.94473442,   12.00100149, ...,   30.62540135,
           30.68166842,   30.73793549],
        [  11.90917952,   11.96544388,   12.02170825, ...,   30.64521234,
           30.7014767 ,   30.75774106],
        [  11.92989169,   11.98615335,   12.042415  , ...,   30.66502333,
           30.72128498,   30.77754664],
        ..., 
        [  30.38443386,   30.43828425,   30.49213464, ...,   48.31661346,
           48.37046385,   48.42431424],
        [  30.40514603,   30.45899372,   30.5128414 , ...,   48.33642445,
           48.39027213,   48.44411982],
        [  30.4258582 ,   30.47970318,   30.53354816, ...,   48.35623544,
           48.41008042,   48.46392539]]])

[ocefpaf]

@kwilcox only @ayan-usgs can answer correctly that question, but you must avoid using rho or psi grids! Otherwise you'll find yourself creating something that is ROMS specific, and that beats the purpose of using pysgrid.

@ayan-usgs I will gladly help you write some docs and/or examples for pysgrid. How far are you from a stable release?

[ayan-usgs]

@kwilcox,

Pysgrid doesn't handle staggered grids as if there are separate grids for nodes for rho and psi points, or for u and v. It instead handles grids that are defined in a variable containing cf_role = grid_topology. For example, in a ROMS dataset, rho, psi, u, v, and etc., would be all on the same grid, defined by an "sgrid" variable. I think it's possible that there could be more than one variable with cf_role = grid_topology within a dataset, but that's not yet supported by pysgrid; I have yet to actually encounter one, so it wasn't very high on priority list.

I was trying to keep method names pretty generic, so .centers returns lon,lat pairs of the grid cell centers, while nodes returns the lon,lat of the nodes.

[rsignell-usgs]

@ayan-usgs, are you using rho, psi, u and v as generic names for the locations on the sgrid, so that these names would apply to locations from any model (e.g. ROMS, WRF, Delft3D...)? That would be okay, but perhaps we could try to see what Arakawa called these locations in his original paper http://www.sciencedirect.com/science/article/pii/B9780124608177500094.

In either case, we should have a way of easily extracting the lon/lat values of the "child" individual curvilinear grids that are associated with the the u, v, psi or rho points, so that it's easy for the user to extract data that will allow them to do analysis and plotting (e.g. pcolormesh(lon_u,lat_u, u))

[kwilcox]

@ocefpaf I agree, I was looking for a generic approach for returning all of the grids. I used the ROMS names for familiarity.

@ayan-usgs @rsignell-usgs Rich said it better, I need a way to access all of the child grids as well as the main grid that we average everything onto for visualization purposes (if necessary). Maybe .centers and .nodes can return the main grid lat/lon pairs, while .grids['some_name'].centers and .grids['some_name'].nodes can return the child grid(s) lat/lon pairs?

[ayan-usgs]

@rsignell-usgs,

No, I am not currently using rho, psi, u, or v as generic terms. Though I've been using centers and nodes as rho and psi, it's the only instance where that occurs; I think of centers and nodes a bit differently as those are the only places where the sgrid object carries data from the original dataset. In general, what currently happens is that the variable name is usually used verbatim. For example, on a ROMS dataset one get get information on that variable by using sg.u, while on a delft3d dataset, the analog would be sg.U1; one could even do sg.lat_rho or sg.lon_rho on a ROMS dataset for instance. However, doing that just gives on attributes of the variable -- none of the data is included.

I think using the original nomenclature that Arakawa used in his paper is an awesome idea. It's something probably be acceptable to many users. It would allow users to work with a dataset without knowing too much about the data; my original thinking was that changing the variable names, when accessing variable attributes, would be confusing to someone who has a dataset and then read it into an sgrid object.

I'll make an issue for this.

[ayan-usgs]

@rsignell-usgs,

So for a child grid, what would be the centers and nodes. For example, in a ROMS dataset would the u variable have lat_u and lon_u be grid centers or grid nodes? My initial thinking is that lat_u and lon_u could be the nodes of the u child grid; these nodes would be coincident with the u vectors. The centers would then be a node of the main grid.

It seems more complicated then it needs to be. Would carrying being able to access lat/lon data via something like sg.u.latitude and sg.u.longitude be sufficient for this need?