lat_ref equals to central_latitude?

[singledoggy]

In the case of wrf, I guess crs should be the crs_new because the mode cent_lat and lon should be in the globe comment here. They are not equal to lat_ref.

@ global String comment MOAD_CEN_LAT =        30.1
@ global String comment STAND_LON =       103

I tested it in one of my WRF out

# central_latitude
central_latitude = 30.1

# new crs that can get the same latlong of wrf
temp = ccrs.LambertConformal(
    globe=globe,
    central_latitude=central_latitude,
    central_longitude=pdef.slon,
    standard_parallels=(pdef.Struelat, pdef.Ntruelat)
)

e, n = pyproj.transform(wgs_proj, temp, pdef.lonref, pdef.latref)

crs_new = ccrs.LambertConformal(
    globe=globe,
    central_latitude=central_latitude,
    central_longitude=pdef.slon,
    standard_parallels=(pdef.Struelat, pdef.Ntruelat),
    false_easting=(pdef.isize - 1) / 2. * pdef.dx - e,
    false_northing=(pdef.jsize - 1) / 2. * pdef.dy - n
)

instead of crs

globe = ccrs.Globe(
    ellipse='sphere', semimajor_axis=6370000, semiminor_axis=6370000)
crs = get_data_projection(ctl, globe=globe)

and the you can get the xcoord and ycoord by lon2d and lat2d via the crs_new, vice versa, Or they may be inconsistent.

lon2d = da.XLONG.values
lat2d = da.XLAT.values

ycoord = np.linspace(0, (pdef.jsize - 1) * pdef.dx, pdef.jsize)
xcoord = np.linspace(0, (pdef.isize - 1) * pdef.dy, pdef.isize)
# case 1
wgs_proj = ccrs.Geodetic(globe=globe)
xx, yy = np.meshgrid(xcoord, ycoord)
x, y, _ = crs.transform_points(wgs_proj, lon2d, lat2d).T
x - xx.T

# %%
# case2
wgs_proj = ccrs.Geodetic(globe=globe)
xx, yy = np.meshgrid(xcoord, ycoord)
x, y, _ = crs_new.transform_points(wgs_proj, lon2d, lat2d).T
x - xx.T

I haven't get a way to automatically get central latitude. And I got the idea from this artical

---

I just use WRF and some other model like NAQPMS, so I'm not sure about the other type of file. If all lat_ref don't equal to central_latitude? Or just WRF need special setting?

[miniufo]

Hi, I am not awared of this point. So I would like to know:

In the case of wrf, I guess crs should be the crs_new because the mode cent_lat and lon should be in the globe comment here. They are not equal to lat_ref.

1. Are crs and crs_new quite different as you initialize them?
2. What is the difference of central_latitude between global string comment and pdef in your WRF case?

[miniufo]

I think basically, WRF and GrADS are two independent softwares. They use somewhat different parameters which would cause inconsistencies (e.g., earth radius is 6371200m in GrADS and 6370000m in WRF; possibly different gravity acceleration g etc.). Maybe WRF can also modify these parameters. I think if you really want a high-level accuracy, you have to construct the crs in exact consistency with WRF, following the pattern provided in utils module. It is a little hard for me to choose one of them as the standard parameters if they are not the same.

This leaves a problem also known as online calculation and offline calculation. Online means you use everything defined in the numerical model while offline means you choose to approximate online in a slightly different way (but you are satisfied with).

Anyway, I would like to parse the global string comment if the WRF ctl provided such information. It is not too hard I guess.

[singledoggy]

1. The difference can be shown below in my case . It's not large, but you can see the edge of the pictures. case1 matches case3.

   
   fig, ax = plt.subplots(1, 3, figsize=(19, 4),
                      subplot_kw={'projection': crs_new})
   da.land.plot(ax=ax[0], transform=crs_new)
   ax[0].title.set_text('crs_new')
   da.land.plot(ax=ax[1], transform=crs)
   ax[1].title.set_text('crs')
   da = da.set_coords(['XLONG', 'XLAT'])
   da.land.plot.pcolormesh(
   ax=ax[2], transform=ccrs.PlateCarree(), x="XLONG", y="XLAT")
   ax[2].title.set_text('LatLon')

<img width="1175" alt="image" src="https://user-images.githubusercontent.com/44222677/166683420-528c2e8a-821b-4988-bc38-4b8798c522f0.png">
3. I can give you my ctl file.

pdef 189 179 lcc 30.100 102.000 95.000 90.000 40.00000 25.00000 103.00000 27000.000 27000.000

and 

@ global String comment MOAD_CEN_LAT = 30.10 @ global String comment STAND_LON = 103.00

https://raw.githubusercontent.com/singledoggy/xgrads/master/ctls/grid.d1.ctl
https://github.com/singledoggy/xgrads/blob/master/ctls/grid.d1?raw=true

[miniufo]

In your case, central_lat (30.1) == MOAD_CEN_LAT (30.1) but lonref (102) != STAND_LON (103). So I guess the bug is at here where pdef.lonref should be changed to pdef.slon.

Also, we may need to add a and b kwargs (6371200) for a perfect sphere in get_coordinates_from_PDEF when construct a Proj in lcc.

Can you make these two small changes and re-plot the figure to see whether it fixes the problem at the edge?

[singledoggy]

central_lat (30.1) == MOAD_CEN_LAT (30.1) is just a coincidence in this case. This doesn't always happen. And I'm pretty sure just change pdef.lonref to pdef.slon won't help. Let's make crs_new2

crs_new2 = ccrs.LambertConformal(
    globe=globe,
    central_latitude=central_latitude,
    central_longitude=pdef.slon,
    standard_parallels=(pdef.Struelat, pdef.Ntruelat),
    false_easting  = pdef.iref * pdef.dx,
    false_northing = pdef.jref * pdef.dy)

I plot (0,0) in each figure to make it clear.

ax[0].plot(0,0, marker="o",markersize=3)
ax[1].plot(0,0, marker="o",markersize=3)
ax[2].plot(0,0, marker="o",markersize=3)

[miniufo]

Can you send me a copy of data (one single time step) and your code snippet for producing the above figure? I will see this in more details. Thank you very much for bringing this issue.

[singledoggy]

Here is my file: https://raw.githubusercontent.com/singledoggy/xgrads/master/ctls/grid.d1.ctl https://github.com/singledoggy/xgrads/blob/master/ctls/grid.d1?raw=true

---

And code: https://raw.githubusercontent.com/singledoggy/xgrads/6e00f304f06a50e116b6e19afa313efa25fa344c/ctls/test.py

[miniufo]

When I run your script, I get an error:

I looked into here and it seems that my pyproj is newer than yours. Do you know how to adapt your code to use the new version of pyproj?

[singledoggy]

I‘m using pyproj==3.3.0 now, and I'm not sure how to adapt it.

[miniufo]

Do you know the difference between CEN_LAT and MOAD_CEN_LAT? Do they differ in certain cases?

[singledoggy]

CEN_LAT and CEN_LON are specific to the nest , and if it's the 2nd or 3rd domain they won't be the same. And lat_ref lon_ref describes the relation of the nest with the Mother of all domains.

[miniufo]

OK, I see. What if the domain is not a nested one? Do we always rely on @ global string comment in ctl file? I guess there are cases that ctl does not contain @ global string comment

[miniufo]

I've fixed this bug in projection. So you can try with the latest codes in get_data_projection. Also, I fixed get_coordinates_from_PDEF accordingly so that the coordinates calculated from this function is very close to the WRF output XLAT and XLONG:

Thank you very much for bringing up this issue.

[singledoggy]

Thanks a lot. It seems that The computation of the domain’s center point is not required for the first WRF domain (it is always centered in its own projection) but it is required for the child domains (which can be placed anywhere in the map).

So if pdef.slon==pdef.lonref and pdef.latref==MOAD_CEN_LAT, then e==n==0. If there is no @ global string comment maybe the users should add it to the ctl files? You may check if pdef.slon==pdef.lonref and raise a warning of the lack information.

[miniufo]

I understand it much better now. So I may need more test cases to remove potential bugs. To me, I think this is deviating from the main purpose of xgrads and would like to stop here. But if you find other issues (maybe using other projections like NPS/SPS), please feel free to report here. Thanks again.

[singledoggy]

Sorry I removed this assignment by mistake. I'm a rookie on Github.