Setting horizontal resolution in pyroms

[sonaljit.m]

Hi all,

I'm trying to setup a ROMS grid for the Puerto Rico and Virgin Islands region using pyroms, and I'm facing an issue with setting the horizontal grid resolution. My objective is to set up a 1km resolution along both zonal and meridional directions, but for some reasons I'm not getting it right. Below is the code that I've written following the Yellow Sea test case (with comments provided):

# Creating the map using Basemap:
# Center point [lat, lon] = [18.2, -65.6]
# Projection is Lambert Conformal,
# with first and 2nd parallels at lat 17.3 and 19.1 respectively.
maproms = Basemap(projection='lcc', lat_0=18.2, lon_0=-65.6, lat_1=17.3, lat_2=19.1,
width=407227, height=199288, resolution='i')
# Zonal and meridional extent are ~400km and ~200km respectively.

#horizontal resolution 1km
delx = 1*1000

# Total number of points=Total distance/grid size (adding 3 extra points as instructed)
Gx = int((maproms.xmax-maproms.xmin)/delx)+3
Gy = int((maproms.ymax-maproms.ymin)/delx)+3

# Setting the domain corners.
lonp = np.array([-67.6, -67.6, -63.73, -63.73])
latp = np.array([19.1, 17.3, 17.3, 19.1])
# Rotation is +1 everywhere. Summation should be 4.
beta = ([1, 1, 1, 1])

# Generate grid
hgrd = pyroms.grid.Gridgen(lonp, latp, beta, (Gx, Gy), proj=maproms)
lonv, latv = maproms(hgrd.x_vert, hgrd.y_vert, inverse=True)
hgrd = pyroms.grid.CGrid_geo(lonv, latv, maproms)

The problem I'm facing is that, the zonal and meridional grid resolution in hgrd are coming out to be different. 1./hgrd.pm shows values around 2000 (2km), whereas 1./hgrd.pn around 500 (0.5km). What I need is to get both 1./hgrd.pm and 1./hgrd.pn to be around 1000.

Probably I'm not entering the variables correctly in gridgen, but I couldn't figure out the error. What am I doing wrong here?

Thanks,
Sonaljit.

[kate]

I don't know why they're so far off, but to do it this way, I would use a Mercator projection. A Lambert conformal conic projection set up this way is not rectangular, but rather a trapezoid (I think) and therefore not orthogonal.

[sonaljit.m]

I tried using Mercator; didn't work. Still the same offset. pm shows ~1/2000 whereas pn ~1/500.

[frederic]

You need to flip Gx and Gy when you call pyroms.grid.Gridgen.
hgrd = pyroms.grid.Gridgen(lonp, latp, beta, (Gy, Gx), proj=maproms)
python use the C matrix indexing and ordering convention.