This is raised by @AminFazlKazemi.

well If the equivalent latitude is chosen perfect, you can fill in values painted as light blue in white space below Equivalent Latitude so The resultant Adiabaticcaly sorted filed is produced.
I chose a mask output from your code corresponding e.g. Latitude=45 degrees North. Pixels are reclassied there with -1 ,+1 and zero. I then multiplied each pixel by its corresponding fractional area in squared kilometers.
If The choosing of Equivalent Latitude is perfect, then integrating area over pixels shown with +1 is exactly equals area over pixels shown with -1. Of course the resolution of dataset plays a role here.

I did not check this carefully. Maybe it is related to the discretized grid points or weights calculated in xgcm. But it needs to be double check.

I made a rough estimate as follows. Grid resolution is 0.7016754 degree. So the area of a grid point at 45°N is (0.7016754*111.11)**2*cos(45°), which is about 4297.32km^2. Similarly, area is 3039.13km^2 at 60N and 1029.35km^2 at 80.25N. It is reasonble to take this as the grid-scale error. I guess your unit of km^2 is wrong (should be m^2 for 10^12).

Also, how do you calculate the area? can you post your code snippet here? Maybe the algorithms are also somewhat different.

Hi

Thanks for your beneficial repo. I just tried to run LWA jupyter notebook. I face these problems:
No module named 'xcontour.xcontour'
I resolved it by substituting from xcontour.xcontour import ... with from xcontour import ...

by the way when I run

dset, grid = add_latlon_metrics(dset)
I face :
NameError: name 'add_latlon_metrics' is not defined

Could you please help

Hi

I am trying to run the LWA code.
I find it difficult install ing geoapps, even thought it is a great package. It seems only we are trying to convert LatLong grid to meters using GeoApps. Is there any problem to the code if we substitute GeoApps function?

for example:

a=6378137.
b=6356752.3
earth_radius=1000.*numpy.sqrt((numpy.power(a*a*numpy.cos(numpy.radians(lat)),2.)+numpy.power(b*b*numpy.sin(numpy.radians(lat)),2.))/(numpy.power(a*numpy.cos(numpy.radians(lat)),2.)+numpy.power(b*numpy.sin(numpy.radians(lat)),2.)))
deg2m=earth_radius*numpy.cos(numpy.radians(lat)) 

After a careful reading of those papers, especially Holliday and McIntyre (1981, JFM), McIntyre and Shepherd (1987, JFM), Roullet and Klein (2009, JFM), Kang and Fringer (2010, JPO), and Huang and Nakamura (2016, JAS), I am pretty sure that local available potential energy (LAPE) is an vertical equivalence to impulse-Casimir wave activity (ICWA) in the horizontal plane. In Huang and Nakamura (2016, JAS), ICWA, and hence LAPE, is clearly demonstrated to be different from local wave activity (LWA). However, in the vertical plane where stratification is usually strong enough to suppress vertical displacement, LWA is quite similar to LAPE in this small-amplitude limit.

So I need to change the calculation of local APE from that of LWA. Local APE is conceptually equivalent to ICWA.

This is raised by @AminFazlKazemi. Let's resolve it one by one.

Any way, I have also another request.
There is also another method for calculating LWA on isentropic surface :
Local Finite-Amplitude Wave Activity as a Diagnostic for Rossby Wave Packets
PAOLO GHINASSI, GEORGIOS FRAGKOULIDIS, AND VOLKMAR WIRTH
Institute for Atmospheric Physics, Johannes Gutenberg University Mainz,
Mainz, Germany
(Manuscript received 26 February 2018, in final form 20 August 2018)
DOI: 10.1175/MWR-D-18-0068.1.

I tried hard to modify you code to this end but in vain.
First : we need to equipartition dm = sigma * ds m being mass .
eqipartitioning of mass not area!
Second we need to integrate (q-Q) sigmads . It is necessary since
Isentropic coordinates will differ peressure sharply (maybe an isentrope is
very low altitude at equator and very high altitude at poles.
This isentropic calculation will add other properties to the resultant LWA.
I would be thankful if you can add this feature to your code.

I see your concern now. You are trying to find the equivalent latitude encompassing the same mass as associated PV contour. I read the Ghinassi's paper. I guess you can try:

1. calculate a discretized mass-PV relation using cal_integral_within_contour, denoted as Mi(Qi);
2. calculate a discretized mass-Lat relation also using cal_integral_within_contour denoted as Mi(Lati) (using latitude contours);
3. Interpolate to remove Mi and get Lati(Qi).

The second point you raised is straighforward when you find Lati(Qi) through Eq. (8) in 2018 paper.

It is interesting to integrate a function that calculates the length of tracer contour and the fractal dimension of it.

• The length of a contour can be calculated using find_contours method in the skimage.measure package, which adopts a marching-cubes algorithm;
• The fractal dimension can be calculated using box-counting algorithm.

This two methods will certainly enrich the xcontour package. Also, the returned results could be a function of equivalent latitudes.