mscross / pysplit

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

Error in trajectory generation on python 3.6.1 #26

Closed zhx215 closed 7 years ago

zhx215 commented 7 years ago

Hi,

I am new to Pysplit, but I used Hysplit before. I read and follow the SCIPY 2015 paper as well as the comments in examples and source code to understand Pysplit technically. When I run the most basic example: bulk_tarajgen_example.py, it shows error message about "can't do nonzero end-relative seeks". Anybody has ideas how to debug this?

In the Command Prompt on Windows 10: I run python bulk_trajgen_example.py, then:

Percent complete: 99.2 Percent complete: 100.0 Complete Hysplit Traceback (most recent call last): File "bulk_trajgen_example.py", line 139, in monthslice=slice(0, 32, 1), get_reverse=True) File "C:\Users\Zhengyu\Anaconda3\lib\site-packages\pysplit\trajectory_generator.py", line 144, in generate_bulktraj meteo_dir, meteofiles, controlfname) File "C:\Users\Zhengyu\Anaconda3\lib\site-packages\pysplit\trajectory_generator.py", line 211, in _reversetraj_whilegen traj.seek(lastline_start, 2) io.UnsupportedOperation: can't do nonzero end-relative seeks

However, I successfully run the bulk_trajgen_example.py on older version of python 2.7 and produce trajectory files. Then I run the basic_plotting_example.py, it got failed as "MatplotlibDeprecationWarning": some of functions in matplotlib, such as get_axis_bgcolor, was not readable in python 2.7. So, I just reinstall the new version of Python to upgrade 3.6.

I need some help to go trough this bugs. Thanks, RZX

mariajmolina commented 7 years ago

I assume other users may be able to provide better help than I can since I am still kind of a python newbie, but I will try in the meantime...

According to the documentation on GitHub, PySPLIT is compatible with Python 2.7 and 3.5... So, I am not sure it would work with Python 3.6. I currently use Python 2.7 with PySPLIT.

Also, after googling your error ("non zero seeks"), the main thing I saw people saying was: "Since the file is a text file, such seeking is not possible" But I have never encountered that error before so I am not sure why that is happening.

It looks like the reason you were getting that "MatplotlibDeprecationWarning" with Python 2.7 was because you were using "get_axisbgcolor", and according to matplotlib that was "Deprecated since version 2.0: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead"._

zhx215 commented 7 years ago

Thanks. I just turned back to Python 2.7, and now I have no problem in running trajectory generation.

Then I try to plot the trajectories by following the example basic_plotting_example.py. Then cmd shows:

C:\Users\Zhengyu\Desktop\pysplit-master\docs\examples>python basic_plotting_example.py C:\Users\Zhengyu\Anaconda2\lib\site-packages\mpl_toolkits\basemap__init.py:1623: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. fill_color = ax.get_axis_bgcolor() C:\Users\Zhengyu\Anaconda2\lib\site-packages\mpl_toolkits\basemap__init.py:1767: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. axisbgc = ax.get_axis_bgcolor() C:\Users\Zhengyu\Anaconda2\lib\site-packages\mpl_toolkits\basemap\init__.py:3260: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0. b = ax.ishold() C:\Users\Zhengyu\Anaconda2\lib\site-packages\mpl_toolkits\basemap\init__.py:3269: MatplotlibDeprecationWarning: axes.hold is deprecated. See the API Changes document (http://matplotlib.org/api/api_changes.html) for more details. ax.hold(b)

Then I edit the basemap init.py: there are two locations using get_axis_bgcolor(), and then I change them to getcolor(). Then I re-run, this error message disappears.

Then I search the init.py codes and find there are 12 locations with ishold function: b = ax.ishold(); how should I edit them? Do I need to change every b = ax.ishold() to ax.hold(b)? Id on't think so (I tried... then NameError: global name 'b' is not defined) Also I did not find any codes with axes.hold. I am not sure what I should do? Anyone has similar experience?

Great thanks.

p.s., I also found if I installed basemap 1.1.0, running basic_plotting_example.py has no outputs, nothing appeared, even no error massages. However, when I degrade to 1.0.7, then error message appears.

mariajmolina commented 7 years ago

I haven't had these issues when plotting and using Python 2.7... Are you modifying the init.py file? I did not have to do that. Can you attach your code as you are running it in python?

mscross commented 7 years ago

Thanks for raising this issue! PySPLIT during trajectory generation does a seek for the last line of a new trajectory file to get initialization information for a reverse trajectory. It looks like the way I used seek is not allowed in Python 3. I will fix how this is done, should be up tomorrow.

Python 3.5 and 2.7 are supported, with 3.6 support planned. I haven't yet tested PySPLIT on 3.6- maybe it works, maybe it doesn't!

Generally, unless you are developing a package, it is not advised to edit its code. If basemap 1.0.7 is throwing errors for you, and 1.1.0 isn't, I suggest using 1.1.0. From your traceback it looks like you have matplotlib 2.0 installed and it doesn't want to play nice with basemap 1.0.7 - the deprecation warning actually refers to matplotlib, not your python version. matplotlib 2.0 appears to be supported for python 2.7.

What environment are you running the example code in? If a Juptyer notebook (recommended), try running %matplotlib inline in one of the cells before generating the figure.

Cheers, Mellissa

zhx215 commented 7 years ago

Hi Mellissa,

Thanks for your guidance. I upgraded to basemap 1.1.0 and now basic_plotting_example.py works on Jupyter notebook. Trajectories could be plotted against its lat and lon. Cool! However, the plot doesn't include the map outline (not drawing the border of countries and coasts). I notice in mapdesigner.py, the class MapDesign defines a function: def init, which includes drawoutlines=True. It means as default the outlines should be plotted as well through basemap.

Do you know what's the reason for this? Here is the code I run in Jupyter notebook:

import pysplit
%matplotlib inline
trajgroup = pysplit.make_trajectorygroup(r'C:/trajectories/colgate/*feb*')
mapcorners =  [-150, 15, -50, 65]
standard_pm = None
bmap_params = pysplit.MapDesign(mapcorners, standard_pm)
bmap = bmap_params.make_basemap()
color_dict = {500.0 : 'blue',
              1000.0 : 'orange',
              1500.0 : 'black'}
for traj in trajgroup:
    altitude0 = traj.data.geometry.apply(lambda p: p.z)[0]
    traj.trajcolor = color_dict[altitude0]
for traj in trajgroup[::4]:
    bmap.plot(*traj.path.xy, c=traj.trajcolor, latlon=True, zorder=20)

image

Great thanks!

mscross commented 7 years ago

I have never seen this before. Its also not drawing the latitude and longitude lines. What happens if you try drawing the countries directly, like bmap.drawcountries()? Also try it with an absurdly high zorder: bmap.drawcountries(zorder=75).

Does this example from basemap work for you?

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
# setup Lambert Conformal basemap.
m = Basemap(width=12000000,height=9000000,projection='lcc',
            resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
# draw coastlines.
m.drawcoastlines()
# draw a boundary around the map, fill the background.
# this background will end up being the ocean color, since
# the continents will be drawn on top.
m.drawmapboundary(fill_color='aqua')
# fill continents, set lake color same as ocean color.
m.fillcontinents(color='coral',lake_color='aqua')
plt.show()
zhx215 commented 7 years ago

Hi Mellissa.

It works. The zorder should be set as high values. Thanks again. Now I could test and produce a variety of plotting. I am still learning the program from example runs. In the example of reversetraj_clippedtraj.py, I just run the same codes provided, then I found the line of warm_trajlist could not create a subset of trajectories with temperature > 0 degree C. When I run codes:

warm_trajlist = [traj for traj in trajgroup if traj.data.Temperature_C[0] > 0.0]
warm_trajgroup = pysplit.TrajectoryGroup(warm_trajlist)
print(trajgroup.trajcount)
print(warm_trajgroup.trajcount)

The jupyter notebook print them as:

360 0

which means there is zero trajectory in the warm_trajlist (every trajectory has start point temperature < 0). I doubt there is something wrong in the along-path environmental variables. How canI call the along-path variables? by just using traj.data.***?

Thanks, RZX

zhx215 commented 7 years ago

oops... I found out. I forgot to check producing along 9 trajectory environmental variables in HYSPLIT after upgrading to 2017 version.

zhx215 commented 7 years ago

I found a bug in the reversetraj_clippedtraj_gen.py:

for traj in warm_trajgroup:
    traj.generate_reversetraj(r'C:/trajectories/umn_example/reversetraj',
                              r'C:/hysplit4/working', r'C:/gdas',
                              meteo_interval='weekly',
                              hysplit="C:\\hysplit4\\exec\\hyts_std")

However, it comes to error message. In traj.py:

def generate_reversetraj(self, hysplit_working, meteo_dir,
                             reverse_dir='default',
                             meteo_interval='weekly',
                             hysplit="C:\\hysplit4\\exec\\hyts_std"):

hysplit_working comes first, meteo comes second, and reverse dir comes third. So the code should be changed to:

for traj in warm_trajgroup:
    traj.generate_reversetraj(r'C:/hysplit4/working',
                              r'C:/gdas',
                              r'C:/trajectories/umn_example/reversetraj',
                              meteo_interval='weekly',
                              hysplit="C:\\hysplit4\\exec\\hyts_std")

Another question I have is when I run traj.generate_clippedtraj(), there are numbers (6 or 7) printed in each traj_generation. I check the source code in traj.py:

if len(meteofiles) == 0:
            raise OSError('No meteorology files found.')
print(len(meteofiles))

I think the number is the len(meteofiles). What's the significant meaning of len(meteofiles) and why is should be printed?

mscross commented 7 years ago

Thanks for bringing that example to my attention! Looks like I forgot to update the example when I changed the code for traj.generate_reversetraj() and traj.generate_clippedtraj() to have a default storage location, rather than requiring the user to provide it.

len(meteofiles) indicates the number of meteorology (GDAS1, in your case) files that PySPLIT found. I probably had it print to serve as a debugging tool. It doesn't need to be printed, so I'll comment that line out.

Yes, you can look at along-path data in the data attribute of a Trajectory. To see what variables you have you can look at the column names, for example for the first Trajectory in your TrajectoryGroup: print(trajgroup[0].data.columns). If one of the column names were Rainfall, then you can print(trajgroup[0].data.Rainfall), etc.

Is your original mapping problem solved?

zhx215 commented 7 years ago

The mapping problem is partly solved; every time I need to add additional two lines at the end:

bmap.drawcoastlines(linewidth=1.0, color='black', zorder=20)
bmap.drawcountries(linewidth=0.5, color='grey', zorder=20)    

Here zorder = 20 is additionally added to make the country and coast outline overlying the white color background; if I set zorder = 15 or lower values, it will be overlaid by white background. I am not sure if this only happens to me. Probably I need to change some parts of code in mapdesigner.py where I see you set by default the outline zorder other low values (~14)?

I think I basically understand the your pysplit program (and python!) and have some ideas to put in my research project. I plan to use it to study the influence of moisture trajectory and moisture uptake source on isotopic composition of precipitation in terrain-complex southernmost South America; it will aid interpretation on downcore isotope data and link isotope with synoptic circulation.