Error running LWA ipynb

[AminFazlKazemi]

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

[miniufo]

Hi, thank you for your interests.

1) the import is a little different on my machine, but it can be fixed as you do; 2) the function is written at here. A quick solution is to copy the function(s) to your code and run with that. This is necessary for xgcm to work properly.

Hope this helps.

[AminFazlKazemi]

Hi I will try and let you be informed ASAP.

On Thu, Apr 28, 2022 at 5:26 AM Yu-Kun Qian @.***> wrote:

Hi, thank you for your interests.

1. the import is a little different on my machine, but it can be fixed as you do;
2. the function is written at here https://github.com/miniufo/GeoApps/blob/bfbac278ff9078d6c2a4ebe6098f7a89e751d7ee/GridUtils.py#L19. A quick solution is to copy the function(s) to your code and run with that. This is necessary for xgcm to work properly.

Hope this helps.

— Reply to this email directly, view it on GitHub https://github.com/miniufo/xcontour/issues/1#issuecomment-1111622484, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMDBA3ICWRO2YJHGTRMB7YLVHHO2BANCNFSM5UQBWVAQ . You are receiving this because you authored the thread.Message ID: @.***>

[miniufo]

I see you raised another issue at #2, but they are the same problem here. You can copy everything you need from GeoApps. Since that package is only a temporary reop, I would not suggest you to install it. Just copy the functions or codes you need into yours will be OK.

[AminFazlKazemi]

Ok, I will do that.

On Thu, Apr 28, 2022 at 12:57 PM Yu-Kun Qian @.***> wrote:

I see you raised another issue at #2 https://github.com/miniufo/xcontour/issues/2, but they are the same problem here. You can copy everything you need from GeoApps. Since that package is only a temporary reop, I would not suggest you to install it. Just copy the functions or codes you need into yours will be OK.

— Reply to this email directly, view it on GitHub https://github.com/miniufo/xcontour/issues/1#issuecomment-1111895526, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMDBA3LJYBUPZDMZDF362ZTVHJDVNANCNFSM5UQBWVAQ . You are receiving this because you authored the thread.Message ID: @.***>

[AminFazlKazemi]

Thanks I managed to run the code using your description.

By the way is there a way I can have Cyclonic and Anti Cyclonic components of LWA too.

Thanks in Advance

On Thu, Apr 28, 2022 at 1:01 PM Amin Fazl Kazemi @.***> wrote:

Ok, I will do that.

On Thu, Apr 28, 2022 at 12:57 PM Yu-Kun Qian @.***> wrote:

I see you raised another issue at #2 https://github.com/miniufo/xcontour/issues/2, but they are the same problem here. You can copy everything you need from GeoApps. Since that package is only a temporary reop, I would not suggest you to install it. Just copy the functions or codes you need into yours will be OK.

— Reply to this email directly, view it on GitHub https://github.com/miniufo/xcontour/issues/1#issuecomment-1111895526, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMDBA3LJYBUPZDMZDF362ZTVHJDVNANCNFSM5UQBWVAQ . You are receiving this because you authored the thread.Message ID: @.***>

[miniufo]

I am not awared that LWA can be decomposed into these two components. Is there any paper talking about this?

[AminFazlKazemi]

Hi equation Number 13 in Huang Nakamura : The Integral over W_ region is Cyclonic LWA The Integral over W+ region is Anti-Cyclonic LWA total LWA is the one you produced in your program. of course it should be divided by latitude cycle length Lx

On Thu, Apr 28, 2022 at 2:17 PM Yu-Kun Qian @.***> wrote:

I am not awared that LWA can be decomposed into these two components. Is there any paper talking about this?

— Reply to this email directly, view it on GitHub https://github.com/miniufo/xcontour/issues/1#issuecomment-1111988857, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMDBA3LJEWSXSB5O7OFXLC3VHJNEBANCNFSM5UQBWVAQ . You are receiving this because you authored the thread.Message ID: @.***>

[miniufo]

Thanks, I'll take a look into it and see if current codes can be refactored slighted to do this. Good idea.

[AminFazlKazemi]

You are welcome. I am the one who is obliged to thank you. I will definitely thank you in my PhD Thesis in Meteorology at Geophysics Institute/University of Tehran/Iran

On Thu, Apr 28, 2022 at 2:41 PM Yu-Kun Qian @.***> wrote:

Thanks, I'll take a look into it and see if current codes can be refactored slighted to do this. Good idea.

— Reply to this email directly, view it on GitHub https://github.com/miniufo/xcontour/issues/1#issuecomment-1112018302, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMDBA3IN4PRBTM3DOMFH4ZDVHJP3JANCNFSM5UQBWVAQ . You are receiving this because you authored the thread.Message ID: @.***>

[miniufo]

Hi, I have added a kwarg part to cal_local_wave_activity, which can be 'all' (default), 'upper' and 'lower', corresponding to all, W+, and W- in Huang and Nakamura (2016, JAS). I feel that upper and lower are somewhat more intuitive than 'cyclonic' and 'anticyclonic', as xcontour also takes into account the vertical cases. I tried with the three parts using their example and verified that all = upper + lower up to machine accuracy. Please give a try and tell me if you find any problem. I am not sure about how the interpretation of the difference seem in the upper and lower parts, but they do look quite different.

By the way, I am also interested in what you are investigating. Do you have any publications that I can get to know?

[AminFazlKazemi]

Hi Many Thanks

My PhD Thesis is about investigating Mediteranean-North Atlantic storm tracks interaction using finite amplitude wave activity diagnostic. There is another simplified form of LWA introduced in this Article :

Wave Events: Climatology, Trends, and Relationship to Northern Hemisphere Winter Blocking and Weather Extremes

PATRICK MARTINEAU AND GANG CHEN Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, Los Angeles, California D. ALEX BURROWS Department of Earth and Atmospheric Sciences, Cornell University, Ithaca, New York (Manuscript received 23 September 2016, in final form 30 January 2017)

DOI: 10.1175/JCLI-D-16-0692.1

Cyclonic& Anticyclonic components both contribute to LWA but every one of them is dominant in some parts of globe as described in this article.

Also It is important to calculate Wave activity 3D flux vectors as formulated in Huang-Nakamura 2016. Formulas Number 20-21-22

About adding a kwarg, I think it is not a good idea. it is needed for the function to create upper and lower LWA components. The user will then simply add them up. In other words : upper_lwa, lower_lwa, contours, masks = analysis.cal_local_wave_activity(tracer, # the original vorticiy ds_latEq.absolute_vorticity, # the sorted vorticity mask_idx=[37,125,170,213]) # select these mask to see the area of LWA

Thanks Again

[miniufo]

There is another simplified form of LWA introduced in this Article

Thank you for your information. I would like to read it.

About adding a kwarg, I think it is not a good idea. it is needed for the function to create upper and lower LWA components. The user will then simply add them up.

I've also considered your suggested way. However, it is clear to me that such way does not speed thing up in current implementation. That is, calculating two components uses almost the same time as calculating two times of a single component. So I just do it in a most simple way by adding a kwarg. If one wants two components, just change the kwarg part for this purpose.

[miniufo]

Also, calculating two components simultaneously would double the memory usage. The current way would be more careful in dealing with the large dataset.

[AminFazlKazemi]

Thanks for your effort. I think calculating Equivalent Latitude is time consuming but maybe it is not so I suggested output separation. Any way, I am going to implement your code and also hnfalwa package of Dr. Huang and test the results. I have also written a code to do this calculation before I became familiar with Dr. Huang package but I put it aside in favor of yours and hers. I'll let you be informed about results comparison between two packages.

[miniufo]

Glad that you are interested in this repo. You could do some practice in your own cases, report some potential bugs, and possibly make contributions.

Actually, Dr. Huang provided many useful help in understanding their theories (see here). Also, I adapted her codes here for my LWA calculation.

Although I am focusing more in oceanic cases, I would like to make this repo more general. So looking forward to your feedbacks.

[AminFazlKazemi]

Sorry for multiple issues. As you know the main superior characteristics of LWA compared to other methods sucjh as Esler/Haynes 1999 lies in equi-partitioning of purturbed PV contour around Lagrangian mean PV corresponding Equivalent Latitude. Well The more we approach to poles and in limiting case, LWA and FAWA approaches zero, By the way we need take into account sphericity of globe in equipartitioning other wise the superior properties of LWA is not guaranteed, namely its validity for finite amplitude waves. """Python

   def earth_radius_function(lats):
       a=6378137.0
       b=6356752.3
       return numpy.sqrt((numpy.power(a*a*numpy.cos(numpy.radians(lats)),2.)+numpy.power(b*b*numpy.sin(numpy.radians(lats)),2.))/(numpy.power(a*numpy.cos(numpy.radians(lats)),2.)+numpy.power(b*numpy.sin(numpy.radians(lats)),2.)))
   earth_radius=earth_radius_function(data.latitude)
   area=(numpy.radians(resolution)*earth_radius)** 2  * (numpy.cos(numpy.radians(data.latitude)))+0 * data.longitude
   for  item_index in numpy.arange(len(masks)):
         item=masks[item_index].copy()
         for field in ['dxC','dyC','rAc']:
             item=item.drop(field)
         item=item.rename({"latitude":"Latitude","longitude":"Longitude"})    
         item=item.assign_coords(latitude=lwa.latitude.values[item_index]).expand_dims("latitude")   
         masks[item_index]=item.copy()
   masks=xr.concat(masks,dim="latitude")    
   ss=masks.sel(latitude=45)
   ss=ss * area
   positive_intrusion_area=xarray.where(ss<0.,ss,0).sum()
   negative_intrusion_area=xarray.where(ss<0.,-ss,0).sum()

"""

The aim of equipartitioning of enclosed area is find an mean lagrangian Potential vorticity level where positive_intrusion_area exactly equals negaive_intrusion_area.

I tried my best in my code to reach to this end. It is more or less true about the outputs of your code but maybe it can be modified.

I suggest or in better words ask you to try it if possible. This will definitely add to the value of your code.

Thanks in advance

[miniufo]

Hi, can you enclose your python codes into ```python codes here ``` for a clear view?

I guess you are talking about the non-sphericity of the earth. But do you know that many numerical models use perfect sphere for earth? Like we ignore the spatial variations of the gravity acceleration g...

[AminFazlKazemi]

Hi

yes, I talk about that but more importantly I talk about exactness of adiabatic PV sorting. I put the code in tripple quotion marks for clarity as you mentioned. 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.

"""python

        vertical_differentiation_method='central'    #'backward'
        if vertical_differentiation_method!='central':
           diffy=-data.pres.diff("level")/data.level.diff("level")/g

 complement=diffy.sel(level=diffy.level.values[0:2]).interp(level=data.level.values[0],method='linear',kwargs={"fill_value":
"extrapolate"})
           diffy=xarray.concat([complement,diffy],'level')

           data['sigma']=diffy.copy().transpose( 'time','level',
'latitude', 'longitude')
        else:

 data['sigma']=(-data['pres'].differentiate("level")/g).transpose(
'time','level', 'latitude', 'longitude')

""" I put here the code to calculate sigma, namely density in isentropic coordinate.

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.

thanks

[miniufo]

Just to be clear,

1. If we assume that the earth is a perfect sphere, the codes in xcontour for estimating area is perfect;
2. xcontour does not provide a function to calculate sigma. One has to do it himself/herself.
3. What is q and Q in your second point? I think the cal_integral_within_contours is enough if you specify integrand=q*Q*sigma as kwargs in isentropic coordinates.

PS: tripple quotion marks you provided for python codes does not take any effect in formating the codes. You can Preview (an bottom at upper left) the results before you submit your comment. Below is a snippet for example (notice the gray background color and keyword highlighted):

import numpy as np

a = np.cos(2)

[AminFazlKazemi]

Well Thanks. I am new to github, learned how to add codes. Thanks I calculated an example : At Latitude=45 degrees North

negative_intrusion_area=1.11084067e+13 squared Kilometers positive_intrusion_area=9.59513011e+12 squared Kilometers

At Latitude=60 degrees North positive_intrusion_area=7.11462961e+12 squared Kilometers negative_intrusion_area=1.55167419e+13 squared Kilometers

At Latitude=80.25 degrees North positive_intrusion_area=2.83549351e+11 squared Kilometers negative_intrusion_area=3.28345966e+12 squared Kilometers

The program works very well of course but may be there exists some way to modify its results.

About my comment, q is potential vorticity and Q is lagrangian mean corresponding Equivalent Latitude I also used your provided kwarg "integrand" as follows : `

analysis = Contour2D(grid, sigma*tracer, dims={'X':'longitude','Y':'latitude'}, dimEq={'Y':'latitude'}, increase=increase,lt=lt)

ctr = analysis.cal_contours(N).rename(var)

mask = xr.where(tracer!=undef, 1, 0).astype(dtype)

table   = analysis.cal_area_eqCoord_table(mask) # A(Yeq) table

area   = analysis.cal_integral_within_contours(ctr,tracer=sigma*tracer,integrand=sigma*0+1).rename('Area')

latEq   = table.lookup_coordinates(area).rename('latEq')

ds_contour = xr.merge([ctr, latEq,area])

preLats = tracer.latitude.astype(dtype)

ds_latEq = analysis.interp_to_dataset(preLats, latEq, ds_contour)

lwa, ctrs, masks = analysis.cal_local_wave_activity(tracer*sigma, ds_latEq[var],mask_idx=numpy.arange(N),part="all")

earth_circle_perimeter=2*numpy.pi*earth_radius*numpy.cos(numpy.radians(lat))

lwa=lwa/earth_circle_perimeter

` But the result is not true . In fact I want to define ctr values based integrand = sigma, also the LatEq and Q be in units of potential vorticity.

but in calculating LWA I want to integrate over integrand=sigma*tracer between contours defined by sigma in first stage.

so finally the calculated LWA and FAWA has the unit "m/s"

I can provide a dataset If you want (sigma and PV_theta).

Again Many Thanks for your efficient code. Sorry for multiple comments

[miniufo]

I calculated an example : At Latitude=45 degrees North negative_intrusion_area=1.11084067e+13 squared Kilometers positive_intrusion_area=9.59513011e+12 squared Kilometers

How do you get these numbers? Which case do you use? and how do you know it is not perfect? (do you have a formula?)

About my comment, q is potential vorticity and Q is lagrangian mean corresponding Equivalent Latitude I also used your provided kwarg "integrand" as follows :

I am a little confused about what you want to do. Could you please write the formula down in a PDF where I can access. I will try to compare the codes above and the formula you want to do and possibly solve your problem.

[AminFazlKazemi]

Well look at the picture at your github about Adiabatic PV sorting at https://github.com/miniufo/xcontour 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 really donot want to bother of course. The only thing I want is to improve your splendid code and make it applicabale to Ghinnassi 2018 case for Isentropic LWA (complement to Nakamura and Solomon 2011 Isentropic FAWA).

The formula is simply Ghinassi 2018 formula 8 for Mass Equivalent Latitude, Equation 15 for Isentropic_LWA.

I am not in a hurry and I think I can close the issue but ask you to do this update in the future to your program.

Thanks

[miniufo]

I open a new issue #3 and continue the discussion.

[Matwayi]

Did you update #GeoApps ? I'm getting this #Rearth error; importError: cannot import name 'Rearth' from 'GeoApps.ConstUtils' (/mnt/g/Azi-Vortex/GeoApps/ConstUtils.py)

[miniufo]

Hi, I've put Rearth into xcontour.utils.py. You can modify your code and import Rearth from there.

[Matwayi]

Fixed! It's ok, thanks for your prompt response.