mscross / pysplit

A package for HYSPLIT air parcel trajectory analysis.
BSD 3-Clause "New" or "Revised" License
151 stars 80 forks source link

moisture uptake and flux analysis #31

Open zhx215 opened 7 years ago

zhx215 commented 7 years ago

Hello, Mellissa

I had some technical questions about moisture uptake and flux parts. I see in the closed issues, Maria has raised it already (Feb 16, 2017), but I was still confused.

moisture = [] for traj in trajgroup: traj.moisture_uptake(precipitation=-0.2,evaporation=0.2,interval=3) moisture.append(traj)

  1. Here she used "moisture.append(traj)" before running "traj.moisture_uptake()". If running "traj.moisture_uptake()", each trajectory will calculate its moisture uptake data (method) and add uptake data to the trajectory (object) itself, is that right? Then, why do I need to use append method to add new traj to moisture list? After running four lines of codes, is the moisture list same with the trajgroup list?

  2. You recommend using scatter plotting to visualize the data with gradient color and you provide example code. After run moisture_uptake(), I can get the information about the uptake location as well as how much moisture is uptaken. Is that correct? Then gradient plotting of moisture uptake is possible. However, as moisture uptake is calculated every 6h, if plotting moisture uptake data along trajectory, do I need to plot against the trajectory midpoint (3h after each window) coordinate within the uptake window?

  3. Any caveat on choosing the evaporation and precipitation threshold value in moisture uptake calculation? I see Baldini et al. (2010) mentioned they used 0.5g/kg specific humidity threshold within 6h. If I want to follow this method, I should set precipitation threshold -0.2, but evaporation threshold 0.7?

  4. For publication, if I want to generate a figure to argue that for example 2010 summer has heavier rainfall isotope because moisture is coming from adjacent ocean while 2011 summer has lighter rainfall isotope because moisture is coming from terrestrial recycled water, I should use plot of moisture uptake or moisture flux? Which one is better in terms of science? I see in your 2015 conference paper, you use moisture flux, but other papers still prefer moisture uptake.

  5. Also, I think in Figure 2, many colored scatter trajectories are overlaid/covered by other trajectories. Have you thought about trying to adding up all trajectory scatter values in one grid and plot something like raster plots? By doing this, you can visualize the total moisture flux/uptake within certain grid in that season, and better say something like the major moisture source is from Bay of Bengal and minor source is from South China Sea. How do you think of this? I see the paper [Extreme rainfall events in southern Sweden: where does the moisture come from?] did this. Another new paper, [Quantification of the impact of moisture source regions on the oxygen isotope composition of precipitation over Eagle Cave, central Spain] plotted moisture uptake percentage on map, which is something I want.

Sorry for too many questions!

RZX

mscross commented 7 years ago

Thanks for the questions!

  1. Putting the trajectories in a list called moisture wasn't strictly necessary. There was no sorting or winnowing of trajectories going on, just a calculation applied to all trajectories. The population of trajectories in moisture is the same as the population of trajectories in trajgroup.

  2. Trajectory.moisture_uptake() produces Trajectory.uptake, which is a GeoDataFramemuch like Trajectory.data: each row comes with a time and location. For Trajectory.uptake, the location is recorded as the midpoint of each window. You can work exclusively with information in Trajectory.uptake to plot uptakes, no need to try to cross reference with Trajectory.data.

  3. Baldini et al's 0.5g/kg threshold is their evaporation threshold. The Sodemann paper with the original algorithm set their evaporation threshold as 0.2 g/kg. A higher evaporation threshold will decrease the sensitivity of the analysis; I believe Baldini et al. discuss the reasoning behind their choice of threshold.

  4. I would choose based on a literature review of similar studies and discussions with your coauthors. My papers on PySPLIT are not scientific studies, but technical articles designed to explain how the package works, in what situations and why is it useful, display its capabilities, and potentially serve as a guide for others seeking to do with other software programs what PySPLIT has done with HYSPLIT.

  5. This can be accomplished with hexbin. In the new version of my paper, recently accepted by Computing in Science and Engineering, I have a figure and some discussion comparing scatter and hexbin plots. A hexbin example is in the works. For what you want, I agree that this type of plot will definitely be more illustrative than scatter plotting a zillion trajectories.

Hope this helps! Mellissa

zhx215 commented 7 years ago

Hi Mellissa,

Thank you. I just took some days to process the moisture uptake data. Now I have a new concern about how to do scatter plotting of the uptake data. For each trajectory, there are some locations with dq_initial high than the evaporation threshold, so these points are identified as the moisture uptake locations. There are two options to plot uptake: one is just to plot the uptake scatter (lon and lat) on the map in the same color; another is to plot scatters in gradient color (the color is determined by dq). Could you give me some guidance on how to do these plotting?

I did gradient color scatter plotting of moisture flux data:

for traj in t_500_rainy_trajgroup:
    mappable = bmap.scatter(*traj.path.xy, c=traj.data.Moisture_Flux.astype(np.float64).values,
                                   cmap=plt.cm.PuBu, vmin=min_moisture_flux,
                                   vmax=max_moisture_flux, zorder=20, latlon=True, edgecolor='none',s = 20)

The difficulty I had is to find the coordinate (x=lon and y=lat) for moisture uptake data (dq_initial > evaporation threshold). In moisture flux data, I can just use *traj.path.xy. How can I do similar thing for traj.uptake? It looks like the xy method is not allowed. 'GeoDataFrame' object has no attribute 'xy'. I am not very into the geopandas and basemap packages. I wish you can give some guidance. Also looking forward to your new moisture uptake example.

mariajmolina commented 7 years ago

Hi Ryan,

When plotting moisture flux data, you are plotting every time step, which coincides with the lat-lon data. Since moisture uptake is computed every 3 or 6 hours (depending what you chose), you need to plot the lat-lon using the same interval. Something like: [::3] or [::6] would do that for your coordinate data. You have to make sure you are coinciding the lat/lon intervals with the moisture uptake data though. Does that make sense?

Maria

zhx215 commented 7 years ago

I found out how to do this: by using geometry & lamda

%matplotlib inline
mapcorners = [-140, -80, -40, -30]
standard_pm = None
bmap_params = pysplit.MapDesign(mapcorners, standard_pm, mapcolor = None, lon_labels=['bottom'], latlon_fs=15)
bmap = bmap_params.make_basemap()

min_dq = 10
max_dq = 0

for traj in t_500_rainy_uptake_trajgroup:
    if traj.uptake.dq.max() > max_dq:
        max_dq = traj.uptake.dq.max()
    if traj.uptake.dq.min() < min_dq:
        min_dq = traj.uptake.dq.min()

for traj in t_500_rainy_uptake_trajgroup:
    mappable = bmap.scatter(traj.uptake.geometry.apply(lambda p: p.x).values,
                            traj.uptake.geometry.apply(lambda p: p.y).values,
                            c=traj.uptake.dq.astype(np.float64).values,
                            cmap=plt.cm.cool,
                            vmin=min_dq,
                            vmax=max_dq,
                            zorder=20,
                            latlon=True,
                            edgecolor='none')
cbar=bmap.colorbar(mappable,ticks=[0.5,1,1.5,2,2.5,3,3.5,4,4.5])
cbar.set_label(r"Moisture Uptake",fontsize=15)

image

mariajmolina commented 7 years ago

It makes way more sense the way you are doing it! I was accessing the regular traj.path coordinates and having to cross reference them, which meant I had to use the intervals.

However, I am now taking a closer look at moisture uptake (I haven't run my experiment yet, I have been just teaching myself how to code) and I am REALLY confused about the traj.uptake.geometry data... When I plot all my traj.uptake dataframe info, the timesteps are correct and going back in time:

                 DateTime  Timestep  Cumulative_Dist  Avg_Pressure  \

Timestep -4.0 2017-01-20 18:00:00 -4.0 NaN 1005.200000 -10.0 2017-01-20 12:00:00 -10.0 NaN 1007.083333 -16.0 2017-01-20 06:00:00 -16.0 NaN 997.700000 -22.0 2017-01-20 00:00:00 -22.0 NaN 979.200000 -28.0 2017-01-19 18:00:00 -28.0 NaN 926.550000 -34.0 2017-01-19 12:00:00 -34.0 NaN 902.666667 -40.0 2017-01-19 06:00:00 -40.0 NaN 814.466667 -46.0 2017-01-19 00:00:00 -46.0 NaN 786.466667 -52.0 2017-01-18 18:00:00 -52.0 NaN 814.333333 -58.0 2017-01-18 12:00:00 -58.0 NaN 777.916667 -64.0 2017-01-18 06:00:00 -64.0 NaN 688.450000 -70.0 2017-01-18 00:00:00 -70.0 NaN 662.816667 -76.0 2017-01-17 18:00:00 -76.0 NaN 652.400000 -82.0 2017-01-17 12:00:00 -82.0 NaN 666.916667 -88.0 2017-01-17 06:00:00 -88.0 NaN 760.216667 -94.0 2017-01-17 00:00:00 -94.0 NaN 779.333333 -100.0 2017-01-16 21:00:00 -100.0 NaN 792.400000

But the traj.uptake.geometry is reversed! Meaning the last lat,lon coordinates are actually the beginning of my backward trajectories. Do you notice this in your output as well?

                                               geometry

Timestep -4.0 POINT Z (-105.767 19.247 2123.6) -10.0 POINT Z (-105.471 19.364 2225.033333333333) -16.0 POINT Z (-105.165 19.658 2251.266666666667) -22.0 POINT Z (-104.443 20.302 2035.016666666666) -28.0 POINT Z (-102.928 21.525 1831.2) -34.0 POINT Z (-101.423 22.555 1553.783333333334) -40.0 POINT Z (-100.366 23.081 1556.483333333333) -46.0 POINT Z (-99.12900000000001 23.525 1439.666666... -52.0 POINT Z (-98.666 23.767 1542.2) -58.0 POINT Z (-98.145 24.113 2039.5) -64.0 POINT Z (-97.67400000000001 24.862 1868.1) -70.0 POINT Z (-96.533 25.461 1005.883333333333) -76.0 POINT Z (-95.124 27.005 759.5166666666668) -82.0 POINT Z (-93.90300000000001 28.541 269.2333333... -88.0 POINT Z (-93.539 29.757 99.08333333333333) -94.0 POINT Z (-93.521 30.288 2.466666666666667) -100.0 POINT Z (-93.20999999999999 30.671 1)

Here is what my geometry looks like in the regular traj.data:

Timestep 0.0 0.0 POINT Z (-92.75 31.56 100) -1.0 0.0 POINT Z (-92.73 31.452 93.3) -2.0 67.3 POINT Z (-92.753 31.334 79.8) -3.0 137.1 POINT Z (-92.821 31.22 34.5) -4.0 208.4 POINT Z (-92.91 31.096 6) -5.0 354.8 POINT Z (-93.00400000000001 30.955 0) -6.0 493.6 POINT Z (-93.099 30.812 0) -7.0 628.1 POINT Z (-93.20999999999999 30.671 0) -8.0 611.4 POINT Z (-93.31 30.549 0) -9.0 593.2 POINT Z (-93.377 30.462 0) -10.0 571.5 POINT Z (-93.423 30.408 0) -11.0 412.0 POINT Z (-93.45999999999999 30.367 0) -12.0 252.1 POINT Z (-93.492 30.327 0) -13.0 92.2 POINT Z (-93.521 30.288 0.7) -14.0 61.4 POINT Z (-93.541 30.245 4.1) -15.0 30.8 POINT Z (-93.554 30.193 10) -16.0 0.0 POINT Z (-93.55500000000001 30.129 22.8) -17.0 0.0 POINT Z (-93.54600000000001 30.04 55.7) -18.0 0.0 POINT Z (-93.541 29.918 94.2) -19.0 0.0 POINT Z (-93.539 29.757 127.2) -20.0 0.0 POINT Z (-93.563 29.571 144.7) -21.0 0.0 POINT Z (-93.611 29.374 149.9) -22.0 0.0 POINT Z (-93.67400000000001 29.167 149.8) ......

mariajmolina commented 7 years ago

Are we supposed to be computing moisture uptake from the end of the backward trajectory?

I need to go back into literature now...

zhx215 commented 7 years ago

Hi Maria, You are right!! The traj.uptake.geometry coordinate data are in reverse to traj.path AND uptake coordinates data. I am analyzing back trajectory. So for me, the traj.uptake.geometry coordinates start from trajectory end point (far/source) to trajectory start point (close/destination); while the uptake data start from "close" to "far". The traj.uptake.geometry was written under so-called windows, which are in reverse to so-called data. Here is the code in traj.py:

windows = self.data.index[::-interval]

So, I need to add [::-1] to make the coordinates correct.

for traj in t_500_rainy_uptake_trajgroup:
        mappable = bmap.scatter(traj.uptake.geometry.apply(lambda p: p.x).values[::-1],
                                traj.uptake.geometry.apply(lambda p: p.y).values[::-1],
                                c=traj.uptake.dq.astype(np.float64).values,
                                cmap=plt.cm.cool,
                                vmin=min_dq,
                                vmax=max_dq,
                                zorder=20,
                                latlon=True,
                                edgecolor='none',
                                s = 100)

Here is the final product: it is much more actually making sense in climatology!! image

Mellissa: I noticed one potential problem in moisture uptake code. In Sodemann et al. (2008a) paper, they did not discount the dq in previous windows after the precipitation h = 0, but in your source code, the dq are also discounted as well by precipitation at h = 0. So I add a line of code to keep the method same with Sodemann et al.:

if self.uptake.loc[w, 'dq_initial'] < precipitation:
                    **if self.uptake.loc[w, 'Timestep'] != 0:**
                        # Adjust previous fractions
                        self.uptake.loc[is_below, 'dq'] = (
                            self.uptake.loc[is_below, 'below'] *
                            self.uptake.loc[w, 'q'])

                        self.uptake.loc[is_above, 'dq'] = (
                            self.uptake.loc[is_above, 'above'] *
                            self.uptake.loc[w, 'q'])

Also, for me, as I want to do scatter plotting of the dq. However, dq also includes above fraction (the moisture source is above boundary layer). I add a line of code in traj.py to exclude the dq above the boundary layer:

        if self.uptake.loc[w, 'dq_initial'] > evaporation:
           ...
            if self.uptake.loc[w, 'above']:
                self.uptake.loc[w, 'dq'] = None

Do you think my editing of code make sense technically?

Thanks, Ryan Z X

mariajmolina commented 7 years ago

Ryan,

I would totally agree with you that we should add a [::-1] to fix the coordinates... but I am not totally convinced because my first coordinate in traj.uptake is (-93.20999999999999 30.671 1) which falls outside of the initial 6 hour interval of my backward trajectory:

Timestep 0.0 0.0 POINT Z (-92.75 31.56 100) -1.0 0.0 POINT Z (-92.73 31.452 93.3) -2.0 67.3 POINT Z (-92.753 31.334 79.8) -3.0 137.1 POINT Z (-92.821 31.22 34.5) -4.0 208.4 POINT Z (-92.91 31.096 6) -5.0 354.8 POINT Z (-93.00400000000001 30.955 0) -6.0 493.6 POINT Z (-93.099 30.812 0) -7.0 628.1 POINT Z (-93.20999999999999 30.671 0)

Shouldn't we be plotting this as the mid-point of every interval?

Maria

mscross commented 7 years ago

The backwards geometry will be fixed in #33 . The midpoint issue is a relic from an older version of pysplit and will also be fixed in that PR. In Trajectory.uptake, there is both a discounted dq (dq) and a non-discounted dq (dq_initial). If you want only the dq_initial below your vertical boundary, there is a way to extract only those values from Trajectory.uptake by taking advantage of the fact that there will be a Nan or something when there isn't an update below your vertical boundary in the is_below or below column. I will post an example. You don't need to edit traj.py to achieve what you want.

Thank you for uncovering these bugs, you guys! I can't believe these slipped through the cracks! I'm currently in the middle of a cross country move and on mobile, but I'll get the geometry issues fixed tonight and work on the example and the 4 digit year thing as I can.

zhx215 commented 7 years ago

Hi Maria,

I think moisture uptake data have the coordinates in the midpoint. In the source code, there is a line of mdpt = interval/2 and [w - mdpt, 'geometry'].x, so mid point locations (and their data/time) are used in geometry.

In my moisture uptake data, the first Timestep is always 0, which corresponds to the geographical location of -3h in trajectory (mid-point; the second is -6, corresponding to -9h., etc.); but in the results you showed, the first one is -4h. I am confused by this. How do you write your the codes for moisture uptake analysis? Copy here and I can have a look to help you.

mariajmolina commented 7 years ago

Ryan,

Here is how I am computing it and getting the -4.0 timestep output ending with -100.0 (the trajectory I am using for testing code has a runtime of -100 hours, and is the one I am using with [0:1:] ).

``

for traj in trajgroup[0:1:]:

traj.moisture_uptake(-0.2,0.5,interval=6,vlim='prs')

print traj.uptake

``

And I get this for geometry (which we already know is backwards):

                                               geometry

Timestep -4.0 POINT Z (-105.767 19.247 2123.6) -10.0 POINT Z (-105.471 19.364 2225.033333333333) -16.0 POINT Z (-105.165 19.658 2251.266666666667) -22.0 POINT Z (-104.443 20.302 2035.016666666666) -28.0 POINT Z (-102.928 21.525 1831.2) -34.0 POINT Z (-101.423 22.555 1553.783333333334) -40.0 POINT Z (-100.366 23.081 1556.483333333333) -46.0 POINT Z (-99.12900000000001 23.525 1439.666666... -52.0 POINT Z (-98.666 23.767 1542.2) -58.0 POINT Z (-98.145 24.113 2039.5) -64.0 POINT Z (-97.67400000000001 24.862 1868.1) -70.0 POINT Z (-96.533 25.461 1005.883333333333) -76.0 POINT Z (-95.124 27.005 759.5166666666668) -82.0 POINT Z (-93.90300000000001 28.541 269.2333333... -88.0 POINT Z (-93.539 29.757 99.08333333333333) -94.0 POINT Z (-93.521 30.288 2.466666666666667) -100.0 POINT Z (-93.20999999999999 30.671 1)

I was comparing my code to yours above and I have a question --> how come you set this:

min_dq = 10 max_dq = 0

Should that be backwards? Like, the max change in q you want is 10? I am currently in a confused state in general.

Thank you for checking!

Maria

mscross commented 7 years ago

Actually, there is no midpoint error. Here's how this algorithm works:

You have a back trajectory; time steps range from 0 to -30. You choose a 6 hour interval. The windowsvariable in the algorithm, which if we printed it would say [-30, -24, -18, -12, -6, 0], indicates at what time step the intervals end.

At -30, the initial conditions are set. From-29 to -24, PySPLIT inspects the conditions during this time and compares it with what happened in the previous interval (in this case just the initial conditions at -30) to determine what has occurred (preciptation, evaporation, indeterminate), and where (above or below the vertical criterion). The midpoint of the interval is time -27, and so PySPLITgoes to Trajectory.data and gets the Point, Cumulative_Dist, and DateTimeat time -27 to represent this interval.

When the -29 to -24 calculation is done, PySPLITinspects -23 to -18, which has a -21 midpoint, performing the same steps as above, except this time adjusting the dq, aboveand belowfraction, and the running totals (unknown, above, below). This repeats, with all previous intervals' dqs, above, below, and totals, adjusted each time. The dq_initialis left alone as a record of the dqprior to adjustment. I will describe the Trajectory.uptake GeoDataFramemore fully in the example.

I hope this is helpful!

mariajmolina commented 7 years ago

THANK YOU SO MUCH! I understand the midpoint thing now and I feel a lot better. I really like this algorithm and look forward to using it.

And just to confirm for comprehension -- the initial conditions are not included in the traj.uptake output, correct? I ask because the first q in traj.uptake is my -4 timestep q from traj.data, which I included images of.

THANKS AGAIN. Much appreciated.

Maria

screenshot 21

screenshot 19

screenshot 20

zhx215 commented 7 years ago

Hi Maria,

Honestly, I don't know why your trajectory starts from -4h. I run the same codes for my HYSPLIT trajectory files.

for traj in t_500_rainy_trajgroup[0:1:]: traj.moisture_uptake(-0.2, 0.5, interval = 6, vlim = 'prs') print (traj.uptake)

                 DateTime  Timestep  Cumulative_Dist  Avg_Pressure  \

Timestep
0.0 2016-12-01 00:00:00 0.0 7.394339e+04 928.066667
-6.0 2016-11-30 18:00:00 -6.0 1.738100e+05 968.683333
-12.0 2016-11-30 12:00:00 -12.0 3.681310e+05 983.700000
-18.0 2016-11-30 06:00:00 -18.0 5.747519e+05 1007.500000
-24.0 2016-11-30 00:00:00 -24.0 7.470974e+05 1011.116667
-30.0 2016-11-29 18:00:00 -30.0 9.173112e+05 1011.816667
-36.0 2016-11-29 12:00:00 -36.0 1.085612e+06 1012.400000
-42.0 2016-11-29 06:00:00 -42.0 1.315419e+06 1013.083333
-48.0 2016-11-29 00:00:00 -48.0 1.600230e+06 1016.066667
-54.0 2016-11-28 18:00:00 -54.0 1.862476e+06 1021.500000
-60.0 2016-11-28 12:00:00 -60.0 2.074800e+06 1024.533333
-66.0 2016-11-28 06:00:00 -66.0 2.236750e+06 1026.600000
-72.0 2016-11-28 00:00:00 -72.0 2.362555e+06 1026.583333
-78.0 2016-11-27 18:00:00 -78.0 2.474321e+06 1027.333333
-84.0 2016-11-27 12:00:00 -84.0 2.587252e+06 1026.933333
-90.0 2016-11-27 06:00:00 -90.0 2.693603e+06 1027.283333
-96.0 2016-11-27 00:00:00 -96.0 2.791191e+06 1027.633333
-102.0 2016-11-26 18:00:00 -102.0 2.892677e+06 1027.733333
-108.0 2016-11-26 12:00:00 -108.0 2.999514e+06 1025.400000
-114.0 2016-11-26 06:00:00 -114.0 3.126810e+06 1024.066667
-120.0 2016-11-26 03:00:00 -120.0 3.191139e+06 1024.000000

I recommend you to re-install the package, keep everything in default, and run some example to see if it is different.

Ryan

mariajmolina commented 7 years ago

Ryan,

Thanks for copying your output so I can compare to mine! I will test this out in a few days. I will let you know if it changes.

Maria

mariajmolina commented 7 years ago

Ryan and Mellissa --

I figured out why my output was starting at -4 timestep!

Its because of the intervals -- When I do the moisture uptake calculation for a 100 hour backward trajectory, the uptake output END starts with -100, but gets cut off at -4. However, when you do the moisture uptake calculation for a 120 hour backward trajectory, the uptake output END starts at -120, and because 120 is divisible by 6, you get to see the 0 timestep. Since 100 is not divisible by 6, it gets cut off at -4. This does not appear to affect the calculation however, it looks like its just a print output thing.

Maria

mariajmolina commented 7 years ago

Now that I think about it further though, only trajectories that are divisible by 6 should probably be used with this algorithm using a 6-hour interval (or divisible by 3 if using a 3-hour interval), otherwise one interval would be weighted differently than the rest...

mscross commented 7 years ago

PySPLIT calculates moisture uptake starting at the end of the trajectory and moving forwards towards time 0. If your trajectory length is not evenly divisible by your interval, it stops before time 0. For example, with the 100 hour trajectory the calculation begins at -100, and moves forward until it has completed all full 6 hour windows, which happens at -4. It does not include or consider time between when it ends and time 0.

However, now with #35 you can choose your starting point (starting_timepoint) for the moisture uptake calculation (Trajectory.moisture_uptake()). PySPLIT will default to the end point of the trajectory. So, with the 100 hr trajectories you can set the starting point to -96, and it will reach time 0. Or you can set it to -30, -77*, whatever you want.

Hope this is helpful!

*This won't reach 0 with a 6 hour interval, of course.

mariajmolina commented 7 years ago

Hi Mellissa,

Thank you so much for taking the time to explain further and adding #35. It is very helpful! I also found the bits below from the Sodemann 2008 paper to be helpful in understanding the algorithm and the backwards/forwards computation parts of it.

[16] 3. Trace the air parcel backward until a positive Δq0 larger than a threshold Δq0c = 0.2 g kg−1 is detected. This threshold suppresses spurious uptakes due to numerical noise and keeps the analysis computationally feasible (see section 5).

[23] 2. Evaluate, proceeding forward in time from the end to start point of the backward trajectory: At an uptake location n inside the BL, calculate the fractional contribution fn of the uptake amount Δqn to the moisture in the air parcel qn as...

Best, Maria

wjjxjd commented 4 years ago

Dear Mellissa,

I think you have done a very good job! Seeing the discussion of the three of you, I learned a lot. But I encountered a problem in visualization, however, there is no problem when testing your code (basic_scatterplotting.py). My code is to imitate the previous one.

`from future import print_function import matplotlib.pyplot as plt import pysplit fig = plt.figure(figsize=(8, 10)) trajgroup = pysplit.make_trajectorygroup(r'/home/wjj/Documents/hysplit1000/20153/*') mapcorners = [-180, -90, 180, -30] standard_pm = [180, -30] param_dict = {'projection':'spstere', 'latlon_labelspacing':(10,30), 'latlon_spacing':(10,15), 'latlon_fs':16, 'resolution':'l', 'mapcolor':'medium','land_alpha':1} bmap_params = pysplit.MapDesign(mapcorners, standard_pm,lon_labels=['bottom'],**param_dict) bmap = bmap_params.make_basemap()

min_dq = 10 max_dq = 0 moisture = [] for traj in trajgroup: traj.moisture_uptake(precipitation=-0.2,evaporation=0.2,interval=6)

moisture.append(traj)

for traj in trajgroup: if traj.uptake.dq.max() > max_dq: max_dq = traj.uptake.dq.max() if traj.uptake.dq.min() < min_dq: min_dq = traj.uptake.dq.min()

for traj in trajgroup: mappable = bmap.scatter(traj.uptake.geometry.apply(lambda p: p.x).values, traj.uptake.geometry.apply(lambda p: p.y).values, c=traj.uptake.dq.astype(np.float64).values, cmap=plt.cm.cool, vmin=min_dq, vmax=max_dq, zorder=20, latlon=True, edgecolor='none') cbar=bmap.colorbar(mappable,ticks=[0.5,1,1.5,2,2.5,3,3.5,4,4.5]) cbar.set_label(r"Moisture Uptake",fontsize=15)`

Traceback (most recent call last): File "mos.py", line 19, in traj.moisture_uptake(precipitation=-0.2,evaporation=0.2,interval=6) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pysplit/traj.py", line 333, in moisture_uptake humidity]]) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 873, in getitem return self._getitem_tuple(key) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 1044, in _getitem_tuple return self._getitem_lowerdim(tup) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 810, in _getitem_lowerdim return getattr(section, self.name)[new_key] File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 879, in getitem return self._getitem_axis(maybe_callable, axis=axis) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 1099, in _getitem_axis return self._getitem_iterable(key, axis=axis) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 1037, in _getitem_iterable keyarr, indexer = self._get_listlike_indexer(key, axis, raise_missing=False) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 1254, in _get_listlike_indexer self._validate_read_indexer(keyarr, indexer, axis, raise_missing=raise_missing) File "/home/wjj/anaconda3/lib/python3.7/site-packages/pandas/core/indexing.py", line 1316, in _validate_read_indexer "Passing list-likes to .loc or [] with any missing labels " KeyError: "Passing list-likes to .loc or [] with any missing labels is no longer supported. The following labels were missing: Index(['Cumulative_Dist'], dtype='object'). See https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#deprecate-loc-reindex-listlike" My python is version 3.7, do you have any suggestions? Thanks!

nitesh1403 commented 3 years ago

Hi, I am getting below errors while calculating moisture uptake. Does anyone have solution for this?

Traceback (most recent call last): File "traj_moist_uptake.py", line 52, in traj.moisture_uptake(precipitation=-0.2,evaporation=0.2,interval=6) File "/opt/miniconda3/envs/pysplitenv/lib/python3.6/site-packages/PySPLIT-0.3.6-py3.6.egg/pysplit/traj.py", line 333, in moisture_uptake humidity]])

. . . .

KeyError: "Passing list-likes to .loc or [] with any missing labels is no longer supported. The following labels were missing: Index(['Cumulative_Dist'], dtype='object'). See https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#deprecate-loc-reindex-listlike"

Thanks in advance.

mscross commented 3 years ago

since you're missing 'Cumulative_Dist' in your trajectory, have you tried running traj.calculate_distance()?

nitesh1403 commented 3 years ago

It worked. Thanks a lot Mellissa. However, still trying to understand the calculations and generated plots. Above discussions helped a lot. I will come again if have more doubts. :)

Thanks again!

-Nitesh

nitesh1403 commented 3 years ago

Hi,

Does anyone calculated total attributed fraction, as Sodemann et. al. 2008?

-Nitesh

mscross commented 3 years ago

Total attributed fraction (i.e., the total attributed to surface uptake) is below_total. 'Below' as a descriptor refers to being under the vlim used to define what is surface uptake and what is not. All of the calculations of Sodemann et al are present, I just chose to name them descriptively- above, below, unknown- rather than use a series of letters (e, f, d).