Merge pull request #14 from xylar/fix_plotting_at_poles

Fix plotting issues by mapping from lon/lat manually

This resolves an unanticipated side-effect of https://github.com/matplotlib/basemap/blob/master/lib/mpl_toolkits/basemap/__init__.py#L510. The coordinates are also being shifted for the option latlon=True before they are converted to map coordinates, viz https://github.com/matplotlib/basemap/blob/master/lib/mpl_toolkits/basemap/__init__.py#L543.  It removes strange plotting errors that occur using standard basemap plotting with the latlon=True option.

The fix produces reasonable plots for standard plot type options, e.g., `['mill2','mill', 'northpole', 'southpole']` with demonstrated success for the following regions:

https://github.com/MPAS-Dev/geometric_features/blob/master/ocean/region/Arctic_Ocean/region.geojson
https://github.com/MPAS-Dev/geometric_features/blob/master/ocean/region/Southern_Ocean/region.geojson
https://github.com/MPAS-Dev/geometric_features/blob/master/landice/region/Antarctica_IMBIE2/region.geojson
https://github.com/MPAS-Dev/geometric_features/blob/master/landice/region/southGreenland1/region.geojson
https://github.com/MPAS-Dev/geometric_features/blob/master/iceshelves/region/Ronne_1/region.geojson

plot_features.py

@@ -107,12 +107,14 @@ def plot_features_file(featurefile, plotname):
 
     fig = plt.figure(figsize=(16,12),dpi=100)
     for anum, maptype in enumerate(['mill2','mill', 'northpole', 'southpole']):
+        print 'plot type: %s'%maptype
         ax = fig.add_subplot(2,2,1+anum)
         feature_num = 0
         for feature in featuredat['features']:
             polytype = feature['geometry']['type']
             coords = feature['geometry']['coordinates']
             feature = feature['properties']['name']
+            print '  feature: %s'%feature
 
             color_num = feature_num % len(colors)
             marker_num = feature_num % len(markers)
@@ -123,20 +125,24 @@ def plot_features_file(featurefile, plotname):
                     for poly in coords:
                         for shape in poly:
                             points = np.asarray(shape)
-                            map.plot(points[:,0], points[:,1], linewidth=2.0, color=colors[color_num], latlon=True)
+                            (x, y) = map(points[:,0], points[:,1])
+                            map.plot(x, y, linewidth=2.0, color=colors[color_num])
                 elif polytype == 'Polygon' or polytype == 'MultiLineString':
                     for poly in coords:
                         points = np.asarray(poly)
-                        map.plot(points[:,0], points[:,1], linewidth=2.0, color=colors[color_num], latlon=True)
+                        (x, y) = map(points[:,0], points[:,1])
+                        map.plot(x, y, linewidth=2.0, color=colors[color_num])
                 elif polytype == 'LineString':
                     points = np.asarray(coords)
                     # due to bug in basemap http://stackoverflow.com/questions/31839047/valueerror-in-python-basemap/32252594#32252594
                     lons = [points[0,0],points[0,0],points[1,0],points[1,0]]
                     lats = [points[0,1],points[0,1],points[1,1],points[1,1]]
-                    map.plot(lons, lats, linewidth=2.0, color=colors[color_num], latlon=True)
+                    (x, y) = map(lons, lats)
+                    map.plot(x, y, linewidth=2.0, color=colors[color_num])
                 elif polytype == 'Point':
                     points = np.asarray(coords)
-                    map.plot(points[0], points[1], markers[marker_num], markersize=20, color=colors[color_num], latlon=True)
+                    (x, y) = map(points[:,0], points[:,1])
+                    map.plot(x, y, markers[marker_num], markersize=20, color=colors[color_num])
                 else:
                     assert False, 'Geometry %s not known.'%(polytype)
             except: