Merge pull request #135 from HHammond/master

Fixed Great Circle Bug, Updated Readme to Markdown.

Author: Jeff Whitaker

Changelog

@@ -4,6 +4,8 @@ version 1.0.8 (not yet released)
   in modification clobbered in upgrade to version 1.1.7 (issue 30).
 * added 'facecolor' keyward argument to drawcounties() method; gives user
   ability to fill counties with specified matplotlib color argument.
+* fix drawgreatcircle bug so that lines exiting and reentering a projection
+  region don't draw horizontally across the map. 
 
 version 1.0.7 (git tag v1.0.7rel)
 ---------------------------------

README README.mdREADME renamed to README.md

@@ -1,28 +1,28 @@
-**Descripton**
+#Basemap
 
-basemap - plot on map projections (with coastlines and political boundaries)
+Plot on map projections (with coastlines and political boundaries)
 using matplotlib.
 
-**Requirements**
+##Requirements
 
-python 2.5 (or higher)
+* python 2.5 (or higher)
 
-matplotlib
+* matplotlib
 
-numpy 
+* numpy 
 
-The GEOS (Geometry Engine - Open Source) library (version 3.1.1 or higher).
+* The GEOS (Geometry Engine - Open Source) library (version 3.1.1 or higher).
 Source code is included in the geos-3.3.3 directory.
 
-PIL (http://pythonware.com/products/pil) is optional (only
+* PIL (http://pythonware.com/products/pil) is optional (only
 needed for Basemap warpimage and bluemarble methods).
 
-On linux, if your python was installed via a package management system, make
+* On linux, if your python was installed via a package management system, make
 sure the corresponding "python-dev" package is also installed.  Otherwise, you
 may not have the python header (Python.h), which is required to build python
 C extensions.
 
-**Copyright**
+##Copyright
 
 source code from proj.4 (http://proj.maptools.org) is included in the
 'src' directory under the terms given in LICENSE_proj4.
@@ -61,7 +61,7 @@ USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
 OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
 PERFORMANCE OF THIS SOFTWARE.
 
-**Documentation** 
+##Documentation 
 
 see http://matplotlib.github.com/basemap/
 
@@ -70,64 +70,42 @@ see scripts in 'examples' directory for example usage.
 read the FAQ and/or email the matplotlib-users mailing list if 
 you have problems or questions.
 
-**Install**
+##Install
 
-0) Install pre-requisite python modules numpy and matplotlib.
+0. Install pre-requisite python modules numpy and matplotlib.
 
-1) Then download basemap-X.Y.Z.tar.gz (approx 100 mb) from
+1. Then download basemap-X.Y.Z.tar.gz (approx 100 mb) from
 the sourceforge download site, unpack and cd to basemap-X.Y.Z.
 
-2) Install the GEOS library.  If you already have it on your
+2. Install the GEOS library.  If you already have it on your
 system, just set the environment variable GEOS_DIR to point to the location 
 of libgeos_c and geos_c.h (if libgeos_c is in /usr/local/lib and
 geos_c.h is in /usr/local/include, set GEOS_DIR to /usr/local).
 Then go to step (3).  If you don't have it, you can build it from
 the source code included with basemap by following these steps:
 
- > cd geos-3.3.3
- > export GEOS_DIR=<where you want the libs and headers to go>
-   A reasonable choice on a Unix-like system is /usr/local, or
-   if you don't have permission to write there, your home directory.
- > ./configure --prefix=$GEOS_DIR 
- > make; make install
+	```
+	 > cd geos-3.3.3
+	 > export GEOS_DIR=<where you want the libs and headers to go>
+	   A reasonable choice on a Unix-like system is /usr/local, or
+	   if you don't have permission to write there, your home directory.
+	 > ./configure --prefix=$GEOS_DIR 
+	 > make; make install
+	```
 
-3) cd back to the top level basemap directory (basemap-X.Y.Z) and
+3. cd back to the top level basemap directory (basemap-X.Y.Z) and
 run the usual 'python setup.py install'.  Check your installation
 by running "from mpl_toolkits.basemap import Basemap" at the python
 prompt.
 
-4) To test, cd to the examples directory and run 'python simpletest.py'.
+4. To test, cd to the examples directory and run 'python simpletest.py'.
 To run all the examples (except those that have extra dependencies
 or require an internet connection), execute 'python run_all.py'.
 
-**Contact**
+##Contact
 
 Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
 
+##Thanks
 
-**Thanks**
-
-to
-
-John Hunter
-Andrew Straw
-Eric Firing
-Rob Hetland
-Scott Sinclair
-Ivan Lima
-Erik Andersen
-Michael Hearne
-Jesper Larsen
-Ryan May
-David Huard
-Mauro Cavalcanti
-Chris Murphy
-Pierre Gerard-Marchant
-Christoph Gohlke
-Eric Bruning
-Stephane Raynaud
-Tom Loredo
-Patrick Marsh
-Phil Elson
-
-for valuable contributions.
+Special thanks to John Hunter, Andrew Straw, Eric Firing, Rob Hetland, Scott Sinclair, Ivan Lima, Erik Andersen, Michael Hearne, Jesper Larsen, Ryan May, David Huard, Mauro Cavalcanti, Chris Murphy, Pierre Gerard-Marchant, Christoph Gohlke, Eric Bruning, Stephane Raynaud, Tom Loredo, Patrick Marsh, Phil Elson, and Henry Hammond for valuable contributions.

examples/nytolondon.py

@@ -17,12 +17,14 @@
 # lonlat, lonlon are lat/lon of London.
 lonlat = 51.53
 lonlon = 0.08
+
 # find 1000 points along the great circle.
 #x,y = m.gcpoints(nylon,nylat,lonlon,lonlat,1000)
 # draw the great circle.
 #m.plot(x,y,linewidth=2)
 # drawgreatcircle performs the previous 2 steps in one call.
 m.drawgreatcircle(nylon,nylat,lonlon,lonlat,linewidth=2,color='b')
+
 m.drawcoastlines()
 m.fillcontinents()
 # draw parallels
@@ -46,6 +48,7 @@
 # lonlat, lonlon are lat/lon of London.
 lonlat = 51.53
 lonlon = 0.08
+
 # find 1000 points along the great circle.
 #x,y = m.gcpoints(nylon,nylat,lonlon,lonlat,1000)
 # draw the great circle.
@@ -62,4 +65,57 @@
 m.drawmeridians(meridians,labels=[1,1,0,1])
 plt.title('Great Circle from New York to London (Gnomonic)')
 sys.stdout.write('plotting Great Circle from New York to London (Gnomonic)\n')
+
+fig=plt.figure()
+
+
+# Example of Great Circles which exit and reenter the map
+m = Basemap(projection='robin', lat_0=0, lon_0=0,resolution='c')
+
+m.drawmapboundary(fill_color='#ffffff',color='#ffffff')
+m.fillcontinents(color='#f2f2f2',lake_color='#ffffff')
+m.drawcoastlines(color='#e0e0e0')
+m.drawcountries(color='#e0e0e0')
+
+parallels = np.arange(-90,90,10.)
+meridians = np.arange(10.,351.,20.)
+m.drawparallels(parallels,labels=[False,False,False,False], color='#d3d3d3',dashes=[1,3])
+m.drawmeridians(meridians,labels=[False,False,False,False], color='#d3d3d3',dashes=[1,3])
+
+# Choose the lat and longtitude of two points
+
+p1 = (45.27,-75.42) # roughly Ottawa
+p2 = (44.05,-4.28)  # roughly Paris
+p3 = (-38.58,145.05) # roughly Victoria
+
+la1, lo1 = p1
+la2, lo2 = p2
+la3, lo3 = p3
+
+# Drawing points; you need to convert from lon-lat to xy-cartesian
+x1,y1 = m(lo1,la1)
+x2,y2 = m(lo2,la2)
+x3,y3 = m(lo3,la3)
+
+# Convert back by setting inverse=True
+lon,lat = m(x1,y1,inverse=True)
+
+# Plot pionts using markers
+m.plot(x1,y1, marker='.', markersize=8, color='#000000')  
+m.plot(x2,y2, marker='.', markersize=8, color='#000000')  
+m.plot(x3,y3, marker='.', markersize=8, color='#000000')  
+
+# Of note, the map uses metres as it's smallest distance, so 1000*x is x km
+# The offset is in projection cooordinates
+plt.text(x1+100000,y1+100000,'Ottawa')
+plt.text(x2+100000,y2+100000,'Paris')
+plt.text(x3+100000,y3+100000,'Victoria')
+
+# Draw a great circle line joining them
+m.drawgreatcircle(lo1,la1,lo2,la2,linewidth=1,color='#000000',alpha=1,del_s=100)
+m.drawgreatcircle(lo2,la2,lo3,la3,linewidth=1,color='#000000',alpha=1,del_s=100)
+
+# Drawing a great circle which exits and reenters the map works on certain projections
+m.drawgreatcircle(lo1,la1,lo3,la3,linewidth=1,color='#FF0000',alpha=1,del_s=100)
+
 plt.show()

lib/mpl_toolkits/basemap/__init__.py

@@ -2891,10 +2891,6 @@ def drawgreatcircle(self,lon1,lat1,lon2,lat2,del_s=100.,**kwargs):
                          method of Basemap instance.
         ==============   =======================================================
 
-        .. note::
-         Cannot handle situations in which the great circle intersects
-         the edge of the map projection domain, and then re-enters the domain.
-
         Returns a matplotlib.lines.Line2D object.
         """
         # use great circle formula for a perfect sphere.
@@ -2908,7 +2904,34 @@ def drawgreatcircle(self,lon1,lat1,lon2,lat2,del_s=100.,**kwargs):
             lats.append(lat)
         lons.append(lon2); lats.append(lat2)
         x, y = self(lons, lats)
-        return self.plot(x,y,**kwargs)
+
+        # Correct wrap around effect of great circles
+
+        # get points
+        p = self.plot(x,y,**kwargs)[0].get_path()
+
+        # since we know the difference between any two points, we can use this to find wrap arounds on the plot
+        max_dist = 1000*del_s*2
+
+        # calculate distances and compare with max allowable distance
+        dists = np.abs(np.diff(p.vertices[:,0]))
+        cuts = np.where( dists > max_dist )[0]
+
+        # if there are any cut points, cut them and begin again at the next point
+        for i,k in enumerate(cuts):
+            # vertex to cut at
+            cut_point = cuts[i]
+
+            # create new vertices with a nan inbetween and set those as the path's vertices
+            verts = np.concatenate(
+                                       [p.vertices[:cut_point, :], 
+                                        [[np.nan, np.nan]], 
+                                        p.vertices[cut_point+1:, :]]
+                                       )
+            p.codes = None
+            p.vertices = verts
+
+        return p
 
     def transform_scalar(self,datin,lons,lats,nx,ny,returnxy=False,checkbounds=False,order=1,masked=False):
         """
@@ -4058,10 +4081,16 @@ def warpimage(self,image="bluemarble",scale=None,**kwargs):
 
         returns a matplotlib.image.AxesImage instance.
         """
+
+        # fix PIL import on some versions of OSX and scipy
         try:
             from PIL import Image
         except ImportError:
-            raise ImportError('warpimage method requires PIL (http://www.pythonware.com/products/pil)')
+            try:
+                import Image      
+            except ImportError:
+                raise ImportError('warpimage method requires PIL (http://www.pythonware.com/products/pil)')
+
         from matplotlib.image import pil_to_array
         if self.celestial:
             msg='warpimage does not work in celestial coordinates'