ESMF regrid problem

[durack1]

From: , Mark
Date: Wednesday, March 4, 2015 at 11:40 AM
To: "Doutriaux, Charles", "Paul J. Durack"
Subject: CDAT regridder issue

Hi Charles and Paul,

Have you guys ever encountered the following issue?

If I regrid using
SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="regrid2”), everything looks okay (see attached top right figure).

But when I regrid using
SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="esmf",regridMethod = "linear”), things get messed up, but only in the NH. Note how things look okay in the SH. (See attached bottom right figure.)

My grids seem to be fine:
Original longitudes:

SWCRE.getLongitude()[:]
Out[39]:
array([ 135.5, 136.5, 137.5, 138.5, 139.5, 140.5, 141.5, 142.5,
143.5, 144.5, 145.5, 146.5, 147.5, 148.5, 149.5, 150.5,
151.5, 152.5, 153.5, 154.5, 155.5, 156.5, 157.5, 158.5,
159.5, 160.5, 161.5, 162.5, 163.5, 164.5, 165.5, 166.5,
167.5, 168.5, 169.5, 170.5, 171.5, 172.5, 173.5, 174.5,
175.5, 176.5, 177.5, 178.5, 179.5, 180.5, 181.5, 182.5,
183.5, 184.5, 185.5, 186.5, 187.5, 188.5, 189.5, 190.5,
191.5, 192.5, 193.5, 194.5, 195.5, 196.5, 197.5, 198.5,
199.5, 200.5, 201.5, 202.5, 203.5, 204.5, 205.5, 206.5,
207.5, 208.5, 209.5, 210.5, 211.5, 212.5, 213.5, 214.5,
215.5, 216.5, 217.5, 218.5, 219.5, 220.5, 221.5, 222.5,
223.5, 224.5, 225.5, 226.5, 227.5, 228.5, 229.5, 230.5,
231.5, 232.5, 233.5, 234.5], dtype=float32)

Destination longitudes:

horiz_grid.getLongitude()[:]
Out[38]:
array([ 135., 137., 139., 141., 143., 145., 147., 149., 151.,
153., 155., 157., 159., 161., 163., 165., 167., 169.,
171., 173., 175., 177., 179., 181., 183., 185., 187.,
189., 191., 193., 195., 197., 199., 201., 203., 205.,
207., 209., 211., 213., 215., 217., 219., 221., 223.,
225., 227., 229., 231., 233., 235.])

Any ideas?

Mark

[doutriaux1]

Can you also send us the latitudes? Or point me to the file.

Thanks

From: "Paul J. Durack" <notifications@github.commailto:notifications@github.com>
Reply-To: UV-CDAT/uvcdat <reply@reply.github.commailto:reply@reply.github.com>
Date: Friday, March 6, 2015 at 4:51 PM
To: UV-CDAT/uvcdat <uvcdat@noreply.github.commailto:uvcdat@noreply.github.com>
Subject: [uvcdat] ESMF regrid problem (#1125)

From: , Mark
Date: Wednesday, March 4, 2015 at 11:40 AM
To: "Doutriaux, Charles", "Paul J. Durack"
Subject: CDAT regridder issue

Hi Charles and Paul,

Have you guys ever encountered the following issue?

If I regrid using
SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="regrid2”), everything looks okay (see attached top right figure).

But when I regrid using
SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="esmf",regridMethod = "linear”), things get messed up, but only in the NH. Note how things look okay in the SH. (See attached bottom right figure.)

My grids seem to be fine:
Original longitudes:

SWCRE.getLongitude()[:]
Out[39]:
array([ 135.5, 136.5, 137.5, 138.5, 139.5, 140.5, 141.5, 142.5,
143.5, 144.5, 145.5, 146.5, 147.5, 148.5, 149.5, 150.5,
151.5, 152.5, 153.5, 154.5, 155.5, 156.5, 157.5, 158.5,
159.5, 160.5, 161.5, 162.5, 163.5, 164.5, 165.5, 166.5,
167.5, 168.5, 169.5, 170.5, 171.5, 172.5, 173.5, 174.5,
175.5, 176.5, 177.5, 178.5, 179.5, 180.5, 181.5, 182.5,
183.5, 184.5, 185.5, 186.5, 187.5, 188.5, 189.5, 190.5,
191.5, 192.5, 193.5, 194.5, 195.5, 196.5, 197.5, 198.5,
199.5, 200.5, 201.5, 202.5, 203.5, 204.5, 205.5, 206.5,
207.5, 208.5, 209.5, 210.5, 211.5, 212.5, 213.5, 214.5,
215.5, 216.5, 217.5, 218.5, 219.5, 220.5, 221.5, 222.5,
223.5, 224.5, 225.5, 226.5, 227.5, 228.5, 229.5, 230.5,
231.5, 232.5, 233.5, 234.5], dtype=float32)

Destination longitudes:

horiz_grid.getLongitude()[:]
Out[38]:
array([ 135., 137., 139., 141., 143., 145., 147., 149., 151.,
153., 155., 157., 159., 161., 163., 165., 167., 169.,
171., 173., 175., 177., 179., 181., 183., 185., 187.,
189., 191., 193., 195., 197., 199., 201., 203., 205.,
207., 209., 211., 213., 215., 217., 219., 221., 223.,
225., 227., 229., 231., 233., 235.])

Any ideas?

Mark
[swcre_regrid_2options]https://cloud.githubusercontent.com/assets/3229632/6539096/fce3ac0c-c420-11e4-8194-e27242995741.png

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/1125.

[mzelinka]

1. Here is what I am using (on feedback):
   $ which cdat
   /usr/local/uvcdat/2015-02-10/bin/cdat
   $ which uvcdat
   uvcdat: aliased to /usr/local/uvcdat/2013-06-05/bin/uvcdat
   $ which ipython
   /usr/local/uvcdat/2015-02-10/bin/ipython

2. Here are the latitudes of the two grids:
   Original latitudes:
   SWCRE.getLatitude()[:]
   Out[3]:
   array([-29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5,
   -20.5, -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5,
   -11.5, -10.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5,
   -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5,
   6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5,
   15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5,
   24.5, 25.5, 26.5, 27.5, 28.5, 29.5], dtype=float32)

Destination latitudes:
horiz_grid.getLatitude()[:]
Out[4]:
array([-29., -27., -25., -23., -21., -19., -17., -15., -13., -11., -9.,
-7., -5., -3., -1., 1., 3., 5., 7., 9., 11., 13.,
15., 17., 19., 21., 23., 25., 27., 29.])

Thanks!

[doutriaux1]

The fact that your uvcdat is different is worrisome (mixed versions) I will try to reproduce next week

[mzelinka]

I just changed to this uvcdat: /usr/local/uvcdat/2015-02-10/bin/uvcdat and the problem remains.

[durack1]

The following code reproduces the error:

import vcs,cdms2
import os,sys
f=cdms2.open("/home/doutriaux1/Mark/charles.nc")
s=f("swcre")

lon_dest = cdms2.createAxis([ 135., 137., 139., 141., 143., 145., 147., 149., 151.,
    153., 155., 157., 159., 161., 163., 165., 167., 169.,
    171., 173., 175., 177., 179., 181., 183., 185., 187.,
    189., 191., 193., 195., 197., 199., 201., 203., 205.,
    207., 209., 211., 213., 215., 217., 219., 221., 223.,
    225., 227., 229., 231., 233., 235.])
lon_dest.designateLongitude()
lon_dest.id="longitude"
lon_dest.units="degrees_east"

lat_dest = cdms2.createAxis([-29., -27., -25., -23., -21., -19., -17., -15., -13., -11., -9.,
    -7., -5., -3., -1., 1., 3., 5., 7., 9., 11., 13.,
    15., 17., 19., 21., 23., 25., 27., 29.])
lat_dest.designateLatitude()
lat_dest.units="degrees_north"
dummy = cdms2.MV2.ones((len(lat_dest),len(lon_dest)))
dummy.setAxisList((lat_dest,lon_dest))

grid_dest=dummy.getGrid()
s.id="orig"
s_regrid2 = s.regrid(grid_dest,regridTool="regrid2")
s_regrid2.id="regrid2"
s_esmf_lin = s.regrid(grid_dest)
s_esmf_lin.id = "ESMF Linear"
s_esmf_con = s.regrid(grid_dest,regridTool="esmf",regridMethod="conservative")
s_esmf_lin.id = "ESMF Conservative"

import EzTemplate

x=vcs.init()
t=x.createtemplate()
t.blank()
t.data.priority=1
t.legend.priority=1
t.dataname.priority=1
t.dataname.y=t.dataname.y*.95
M=EzTemplate.Multi(template=t,x=x,rows=2,columns=2)
gm=x.createboxfill()
levels= vcs.mkscale(*vcs.minmax(s))
cols = vcs.getcolors(levels)
gm.boxfill_type = "custom"
gm.fillareacolor = cols
gm.levels = levels
x.plot(s,M.get(),gm)
x.plot(s_regrid2,M.get(),gm)
x.plot(s_esmf_lin,M.get(),gm)
x.plot(s_esmf_con,M.get(),gm)
x.png("esmf_issue_1125")
raw_input("Press enter")

And here is the output:

[durack1]

@dnadeau4 @doutriaux1 just pinging you both on this..

[doutriaux1]

it's on my list for 2.4... keep your finger crossed....

[doutriaux1]

@dnadeau4 this is the answer I got from Ryan OKuinghttons

Hi Charles,

I've got some good news and some bad news for you.  The good news is 
that I've found the issue, and the bad news is that it was something in 
my ESMP reproducer script.  This means that your problem is likely 
caused by the wrapper layer between ESMP and UVCDAT.  The mistake I made 
in the reproducer was caused by creating a periodic Grid, instead of a 
regional Grid.  This is easy to fix by replacing:

grid = ESMP.ESMP_GridCreate1PeriDim(maxIndex)

with:

grid = ESMP.ESMP_GridCreateNoPeriDim(maxIndex)

After this change, and a few other minor mods to clean up the plots and 
make the two reproducer codes more alike, the results are as expected.  
I attached the ESMP and ESMPy reproducer codes, you can easily reproduce 
this error by reversing the change I mentioned above.  There is also one 
minor error in the ESMPy script that I still need to fix, but I wanted 
to get you a more satisfactory solution as soon as possible.

If I can make a guess, I would say that the error in the interface layer 
(as well as the one I made in the reproducer) is due to the fact that 
the latitude data in this example _appears_ to cross a pole (at 0 
degrees).  I suppose that the grid in this example is intended to have 
the poles at -90 and 90 degrees.  Either way, I think that the structure 
of this data is throwing off the interface layer, causing it to create a 
Grid with a periodic dimension. Again, this is pure speculation on my part.

Also, when talking to Ben today I realized that I forgot to mention in 
my last mail that ESMPy now also interfaces with OCGIS using converters 
between ESMPy Fields and OCGIS Variables.  The current capability is 
limited to IO and subsetting in serial, but we plan to expand that in 
the future to include computations and parallel execution.

Thanks!

Ryan

[doutriaux1]

Ryan scripts:
uvcdat_error_esmp.py

# This is an ESMP reproducer of uvcdat bug 1125, esmf support request 3613723
import ESMP
import numpy as np

def make_index( nx, ny ):
    return np.arange( nx*ny ).reshape( ny, nx )

def create_grid(lons, lats, lonbnds, latbnds):

    maxIndex = np.array([len(lons), len(lats)], dtype=np.int32)

    grid = ESMP.ESMP_GridCreateNoPeriDim(maxIndex)

    ##   CENTERS
    ESMP.ESMP_GridAddCoord(grid, staggerloc=ESMP.ESMP_STAGGERLOC_CENTER)

    exLB_center, exUB_center = ESMP.ESMP_GridGetCoord(grid,
                                                      ESMP.ESMP_STAGGERLOC_CENTER)

    # get the coordinate pointers and set the coordinates
    [x, y] = [0, 1]
    gridXCenter = ESMP.ESMP_GridGetCoordPtr(grid, x, ESMP.ESMP_STAGGERLOC_CENTER)
    gridYCenter = ESMP.ESMP_GridGetCoordPtr(grid, y, ESMP.ESMP_STAGGERLOC_CENTER)

    # make an array that holds indices from lower_bounds to upper_bounds
    bnd2indX = np.arange(exLB_center[x], exUB_center[x], 1)
    bnd2indY = np.arange(exLB_center[y], exUB_center[y], 1)

    pts = make_index(exUB_center[x] - exLB_center[x],
                     exUB_center[y] - exLB_center[y])

    for i0 in range(len(bnd2indX)):
        gridXCenter[pts[:, i0]] = lons[i0]

    for i1 in range(len(bnd2indY)):
        gridYCenter[pts[i1, :]] = lats[i1]

    ##   CORNERS
    ESMP.ESMP_GridAddCoord(grid, staggerloc=ESMP.ESMP_STAGGERLOC_CORNER)

    exLB_corner, exUB_corner = ESMP.ESMP_GridGetCoord(grid, ESMP.ESMP_STAGGERLOC_CORNER)

    # get the coordinate pointers and set the coordinates
    [x,y] = [0, 1]
    gridXCorner = ESMP.ESMP_GridGetCoordPtr(grid, x, ESMP.ESMP_STAGGERLOC_CORNER)
    gridYCorner = ESMP.ESMP_GridGetCoordPtr(grid, y, ESMP.ESMP_STAGGERLOC_CORNER)

    # make an array that holds indices from lower_bounds to upper_bounds
    bnd2indX = np.arange(exLB_corner[x], exUB_corner[x], 1)
    bnd2indY = np.arange(exLB_corner[y], exUB_corner[y], 1)

    pts = make_index(exUB_corner[x] - exLB_corner[x],
                     exUB_corner[y] - exLB_corner[y])

    for i0 in range(len(bnd2indX)-1):
        gridXCorner[pts[:, i0]] = lonbnds[i0, 0]
    gridXCorner[pts[:, i0+1]] = lonbnds[i0, 1]

    for i1 in range(len(bnd2indY)-1):
        gridYCorner[pts[i1, :]] = latbnds[i1, 0]
    gridYCorner[pts[i1+1, :]] = latbnds[i1, 1]

    return grid

def build_analyticfield(field, grid):
    # get the field pointer
    fieldPtr = ESMP.ESMP_FieldGetPtr(field)

    # get the grid bounds and coordinate pointers
    exLB, exUB = ESMP.ESMP_GridGetCoord(grid, ESMP.ESMP_STAGGERLOC_CENTER)

    # make an array that holds indices from lower_bounds to upper_bounds
    [x,y] = [0, 1]
    bnd2indX = np.arange(exLB[x], exUB[x], 1)
    bnd2indY = np.arange(exLB[y], exUB[y], 1)

    realdata = False
    try:
        import netCDF4 as nc
        f = nc.Dataset('charles.nc')
        swcre = f.variables['swcre']
        swcre = swcre[:]
        realdata = True
    except:
        raise ImportError('netCDF4 not available on this machine')

    p = 0
    for i1 in range(len(bnd2indX)):
        for i0 in range(len(bnd2indY)):
            if realdata:
                fieldPtr[p] = swcre.flatten()[p]
            else:
                fieldPtr[p] = 2.0
            p = p + 1

    return field

def compute_mass(valuefield, areafield, fracfield, dofrac):
    mass = 0.0
    ESMP.ESMP_FieldRegridGetArea(areafield)
    area = ESMP.ESMP_FieldGetPtr(areafield)
    value = ESMP.ESMP_FieldGetPtr(valuefield)
    frac = 0
    if dofrac:
        frac = ESMP.ESMP_FieldGetPtr(fracfield)
    for i in range(valuefield.size):
        if dofrac:
            mass += area[i]*value[i]*frac[i]
        else:
            mass += area[i]*value[i]

    return mass


def plot(srclons, srclats, srcfield, dstlons, dstlats, interpfield):

    try:
        import matplotlib
        import matplotlib.pyplot as plt
    except:
        raise ImportError("matplotlib is not available on this machine")

    exact = ESMP.ESMP_FieldGetPtr(srcfield)
    exact = exact.reshape(len(srclats.flatten()), len(srclons.flatten()))

    interp = ESMP.ESMP_FieldGetPtr(interpfield)
    interp = interp.reshape(len(dstlats.flatten()), len(dstlons.flatten()))

    fig = plt.figure(1, (15, 6))
    fig.suptitle('ESMP Conservative Regridding', fontsize=14, fontweight='bold')

    ax = fig.add_subplot(1, 2, 1)
    im = ax.imshow(exact, vmin=-140, vmax=0, cmap='gist_ncar', aspect='auto',
                   extent=[min(srclons), max(srclons), min(srclats), max(srclats)])
    ax.set_xbound(lower=min(srclons), upper=max(srclons))
    ax.set_ybound(lower=min(srclats), upper=max(srclats))
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_title("Source Data")

    ax = fig.add_subplot(1, 2, 2)
    im = ax.imshow(interp, vmin=-140, vmax=0, cmap='gist_ncar', aspect='auto',
                   extent=[min(dstlons), max(dstlons), min(dstlats), max(dstlats)])
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_title("Conservative Regrid Solution")