from mpl_toolkits.basemap import Basemap
m = Basemap()
m.bluemarble()
m.drawcoastlines()

The blue marble map is upside down, both on screen and in relation to the coastlines. Worked before :)

matplotlib and basemap from git (2012-06-25). PIL 1.1.7

Mac osx 10.7, with python through homebrew.

The script below runs just fine with basemap version 1.0.1, but crashes with the following error on version 1.0.2.:
Traceback (most recent call last):
File "./maptest.py", line 33, in
linewidth=1,color='black',yoffset=yoff,xoffset=xoff,dashes=[1,0])
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/mpl_toolkits/basemap/init.py", line 2067, in drawparallels
if t is not None: linecolls[int(lat)][1].append(t)
KeyError: 38

This script is a trivialized snippet of a real function, where the meridian and parallel values are dynamically derived. The instance that crashed generated the values below.

Script:

!/usr/bin/env python

import matplotlib

use the non-interactive matplotlib setting

matplotlib.use('agg')
from mpl_toolkits.basemap import Basemap
from pylab import *

figheight = 5.35891
figwidth = 5.4
bounds = [-89.704166666702776, -85.962500000037608, 36.987499999949222, 39.895833333281395]
clat = bounds[2] + (bounds[3] - bounds[2])/2
clon = bounds[0] + (bounds[1] - bounds[0])/2
fig = figure(figsize=(figwidth,figheight),edgecolor='g',facecolor='g')
ax1 = fig.add_axes([0,0,1.0,1.0])
mapcontour = Basemap(llcrnrlon=bounds[0],llcrnrlat=bounds[2],
urcrnrlon=bounds[1],urcrnrlat=bounds[3],
resolution='h',projection='merc',lat_ts=clat)

par = array([ 37. , 37.75, 38.5 , 39.25])
mer = array([-89.5, -88.5, -87.5, -86.5])
xmap_range = mapcontour.xmax-mapcontour.xmin
ymap_range = mapcontour.ymax-mapcontour.ymin
xoff = -0.09_(xmap_range)
yoff = -0.04_(ymap_range)
mapcontour.drawmeridians(mer,labels=[0,0,1,0],fontsize=8,
linewidth=1,color='black',yoffset=yoff,xoffset=xoff,dashes=[1,0])
mapcontour.drawparallels(par,labels=[0,1,0,0],fontsize=8,
linewidth=1,color='black',yoffset=yoff,xoffset=xoff,dashes=[1,0])
mapcontour.drawmapboundary(color='k',linewidth=2.0)
savefig('test.png')

I've been having a problem in my basemap plots of areas with small spatial extents, where my meridians are cutoff like so:

I've tried simplifying the code as much as possible, but still get a similar problem:

The meridians plot correctly when I use a normal mercator instead of the transverse mercator projection.

Here's my simplified code:

import numpy as np
import matplotlib as mpl
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

"""Input desired lats (y) and lons(x) of plot"""

xmin = -166.22
xmax = -166.134

ymin = 53.7439
ymax = 53.7827

m = Basemap(projection='tmerc', lon_0=-165. , lat_0=0 , llcrnrlon=xmin, llcrnrlat=ymin, urcrnrlon=xmax , urcrnrlat=ymax)

meridians = np.arange(xmin,xmax,(30/3600)) # 30 arcsecond spacing
m.drawmeridians(meridians)

parallels = np.arange(ymin,ymax,(30/3600))
m.drawparallels(parallels)

plt.show()

Hi:

I am not sure if it worth your attention. But I am developing in pycharm IDE and previously everything worked fine with basemap. Just recently I saw that it is not able to find the installation and I am not able to use autocompletion and other conveniences.

I have created an issue on pycharm website:

http://youtrack.jetbrains.com/issue/PY-13262

If you think there is something that could be done from basemap side to facilitate integration it would be really cool.

Thanks

I am a little bit hesitant about the name "arcgisimage" for the new feature added to basemap recently. First off, I am concerned about using a trademark in a function name. Second, is it over-specifying what it is that this function can do? Why not "wmsimage"?

Lastly, on a related note, does basemap have an rcparams mechanism for storing defaults like the url and service names?

Hello,

It appears that the change for shapefile.py and the default_encoding being able to be set got clobbered with a merge. That means that drawcounties is broken again in the latest release. I'm fixing mine locally but figured I would let you know.

It appears that drawgreatcircle() has some sign issues, especially when an equator crossing is included.

I've attached a picture showing the resulting paths created by drawgreatcircle() using a few major cities. I called great circle in both directions for each pair of city. It appears that it only works if the direction is west-to-east and no equator crossing is involved. The Tokyo-to-LA path shows the expected behavior, so crossing 180 longitude is ok.

Ignore the red to in the middle of the Atlantic Ocean, which I put there manually for testing.

My basemap version is whatever is in the Ubuntu 12.10 standard repositories.

I am trying to scatter plot on maps with nonzero middle longitude, but the areas where the map is "wrapped" are not plotted.

My test script:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

m = Basemap(lon_0=160)
m.fillcontinents()

for x in range(-179, 179):
px1, py1 = m(x, x/4.0)
m.scatter(px1,py1, zorder=3)

plt.show()

If lon_0 is zero it works fine but with non-zero lat_0 some areas are not plotted.

It seems like the m(lon,lat) is not wrapping the points, the points I expect to wrap are negative in px1?

Or am I missing some transformation step?

When using an equidistant cylindrical projection over a selected subregion of the globe, e.g. map = Basemap(llcrnrlon=50,llcrnrlat=-80,urcrnrlon=300,urcrnrlat=70,projection='cyl'), any contour plot (i.e. any plot generated using the map.contour() function) cannot plot any values south of 10N. In other words, everything south of 10N just shows up blank on the plot.

Meridians and parallels don't display properly on a map of a very
small area because lats and lons arrays are generated with a fixed
step of 0.01.
Is it possible to adapt this step to the current map (for instance
using attributes like urcrnrx, etc, as with xoffset and yoffset)?

Using the following setup, Python 2.7/matplotlib 1.2.1/basemap 1.0.6

from mpl_toolkits import basemap

m1 = basemap.Basemap(projection='cea', lon_0=0)
lon = [-135, -45, 45, 135]
lat = [45, 45, 45, 45]

If I then run the following command, the line plots as expected:

m1.plot(lon, lat, latlon=True)

But plotting just the first two points alone gives a failure:

m1.plot(lon[:2], lat[:2], latlon=True)

Here's the result:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-273-c6f75eb64ab8> in <module>()
      8 lat = [45, 45, 45, 45]
      9 #m1.plot(lon, lat, latlon=True)
---> 10 m1.plot(lon[:2], lat[:2], latlon=True)

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in with_transform(self, x, y, *args, **kwargs)
    525             if self.projection in _cylproj or self.projection in _pseudocyl:
    526                 if x.ndim == 1:
--> 527                     x = self.shiftdata(x)
    528                 elif x.ndim == 0:
    529                     if x > 180:

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in shiftdata(self, lonsin, datain, lon_0)
   4700             londiff = np.abs(lonsin[0:-1]-lonsin[1:])
   4701             londiff_sort = np.sort(londiff)
-> 4702             thresh = 360.-londiff_sort[-2]
   4703             itemindex = len(lonsin)-np.where(londiff>=thresh)[0]
   4704             if itemindex:

IndexError: index -2 is out of bounds for axis 0 with size 1

Edit: here's an error from a similar call, but I believe it's unrelated to the first:

lat = [-45, 45, 45, 45, 45, 0]
lon = [135, -135, -45, 45, 135, -90]
m1.plot(lon, lat, latlon=True)

error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-295-d815d338a0b6> in <module>()
     11 lon = [135, -135, -45, 45, 135, -90]
     12 
---> 13 m1.plot(lon, lat, latlon=True)
     14 

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in with_transform(self, x, y, *args, **kwargs)
    525             if self.projection in _cylproj or self.projection in _pseudocyl:
    526                 if x.ndim == 1:
--> 527                     x = self.shiftdata(x)
    528                 elif x.ndim == 0:
    529                     if x > 180:

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in shiftdata(self, lonsin, datain, lon_0)
   4702             thresh = 360.-londiff_sort[-2]
   4703             itemindex = len(lonsin)-np.where(londiff>=thresh)[0]
-> 4704             if itemindex:
   4705                 # check to see if cyclic (wraparound) point included
   4706                 # if so, remove it.

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I'm trying to zoom in on a cylindrical projection of some data and think there's a Basemap 1.0.5 issue. The data I'm plotting is from a climate model and is of fairly coarse resolution. The plot at http://climate.ncas.ac.uk/ajh/basemap/python.png shows the problem with missing contours although there is data there after applying addcyclic and shiftgrid.

The data in longitude starts at 0 degrees and then has steps of 3.75 degrees to 356.25 degrees. When the plotting region starts or ends at a point that isn't a grid point location in longitude the plot shows this issue. If I change the contour region to -7.5 to 3.75 (both of which are grid points) then the contours are all drawn correctly.

The simple code and data is available at http://climate.ncas.ac.uk/ajh/basemap. is this a known issue?

Regards
Andy

Hi ladies and gentlemen,

I was wondering how the *.dat files were generated, and if it would be possible to annotate some of them with their names in some standardized format. For instance, for countries, we could append another field after the #bytes value that contains the ISO2 name of the country of the coordinates listed. Then one could build LineCollections for each country and do fun per-country graphing.

Otherwise, we need to use shapefiles -- not all that bad except they might not match the specified resolution or even the data included in basemap.

I'd love to submit patches, but unfortunately there's no rhyme or reason to the internal organization of the *.dat files. I hypothesize having their source will facilitate such a transformation.

Thoughts?

In drawmapboundary the line/edge color of the limb is set in all cases except in the round=True case. This should not be.

The following is probably the problematic line:

            limb = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
                    radius=0.5*(self.xmax-self.xmin),fc='none')

Hi,

I am not sure if this is the right place to ask, so feel free to send me elsewhere with the question. I have created a map plot and a colorbar. For the colorbar I've used the format argument, supplying the following object as a value:

#format the colorbar tick labels
sfmt = ScalarFormatter(useMathText=True)
sfmt.set_powerlimits((-2, 2))

As you can see on the image below my text from the multiplier is overlapping with a top tick label. I would like to lift it a bit, is there an easy way?
I tried this, but it does not work, since after basemap.colorbar() the title does not contain the multiplier yet:

ax = colorbar.ax
title = ax.get_title()
ax.set_title("{0}\n\n\n\n---".format(title))

I could, probably increase the width of the colorbar or decrease font size... But I like this width and the font size seems to be very readable.

Thanks

Sasha

Hi

I am making Polar Stereographic Projection maps of some climate model outputs. For some of these data, the plot looks weird. For example, in this figure
(https://dl.dropboxusercontent.com/u/45427012/ALB_ERR.png), only two color contours showed up while the actual data should span much wider range. Furthermore, a large portion of the region should be blank since the data are masked out by netcdf module already (they are undefined).
You can find the script and dataset which generated the figure at here

python script
https://dl.dropboxusercontent.com/u/45427012/plot_albice_error.py

data (netcdf4), you'll need netcdf module to open it.
https://dl.dropboxusercontent.com/u/45427012/alb.nc4

BTW, I am using basemap-1.0.6 with matplotlib-1.2.1 on py2.7.

Any help will be much appreciated.

Bin Zhao

For the labelling in drawmeridians/drawparallels the docs state:

labels    list of 4 values (default [0,0,0,0]) that control whether meridians
          are labelled where they intersect the left, right, top or bottom of
          the plot. For example labels=[1,0,0,1] will cause meridians to be
          labelled where they intersect the left and and bottom of the plot,
          but not the right and top.

However for the following example this doesnt seem to work as expected (with current master):

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

bmap = Basemap(projection='aeqd', resolution="i", area_thresh=1000.0,
               lat_0=39.9, lon_0=58.7, width=4378000, height=1223000)

meridians = [35.,  40.,  45.,  50.,  55.,  60.,  65.,  70.,  75.,  80.]
bmap.drawmeridians(meridians, labels=[1, 0, 0, 1])
parallels = [35.,  40.,  45.]
bmap.drawparallels(parallels, labels=[0, 1, 1, 0])

plt.show()

Note the parallels labels at the top and the one meridian label on the left side.

Hi

I'm trying to plot vectors in the Southern ocean using basemap and matplotlib.
To illustrate my problem I've created a zonal jet at 60S.
I have U and V as 2D arrays and lon and lat as 1D arrays, all with cyclic point added and lon is from -180 to 180.
If I plot in orthographic projection the vectors correctly shows eastward motion, but if I use any of the polar projections (spstere, splaea spaeqd) the vectors reverse!
This is really strange and oddly enough it does not seem to happen if I make the same plot over the North Pole.

I've noted that in the southern hemisphere lon_0 has a different meaning for orthographic projection than for spstere, splaea etc.
So to orient the maps the same way I must pass lon_0=0 in orthographic projection and lon_0=180 in polar stereographic projection if I plot southern hemisphere.
But this is not the case in the Northern hemisphere.
Could this be part of the problem?

Here's the script that creates the results:
https://dl.dropboxusercontent.com/u/86754720/plot_quiver_bug_test.py
Here's the contourf of u and v
https://dl.dropboxusercontent.com/u/86754720/velocity_bug_contourf.png
and here are the vectors
https://dl.dropboxusercontent.com/u/86754720/velocity_bug_quiver.png

Hello,

First of all, thank you for the great module! I find it incredibly useful.

I surveyed previous issues and haven't found one corresponding to my problem. Here is my issue:
Using Basemap, my colleagues and I noticed that some islands tend to be missing from the map. Given the amount of detail in the mainland shore nearby, it doesn't appear to be a resolution problem. They seem to be missing when using drawcoastlines and fillcontinents but they do appear (although without coastlines) when using drawlsmask. The first missing island we noticed is the (entire) Canadian province of Prince Edward Island. Then, we also noticed the absence of Anticosti Island, and others. They are fairly big (see figures).

Is that intended? Let me know if I can be of any help investigating this.

Thanks!
Steph

Eric Firing reported on mpl-developer:

"If I do a clean install of mpl master, and then of basemap, basemap lands in dist-packages/mpl_toolkits, as it always has. But now it is not found--I can't import it. It seems that now the real mpl_toolkits is cleverly hidden inside an egg directory with a monstrous name, leaving basemap orphaned in a directory with no init.py. As a workaround I can symlink it into the egg location. "

I can see three potential solutions for this problem:

1. just add an empty init__py to basemap in lib/mpl_toolkits. Advantage: easy, and allows user code with "from mpl_toolkits.basemap import Basemap" to still work. Disadvantage: packaging systems will have to deal with the potential of colliding mpl_toolkits/__init.py files (one owned by matplotlib and one owned by basemap).

2. somehow figure out how to install inside a matplotlib egg, or to make a symlink to the mpl_toolkits/init.py inside the matplotlib egg. This is probably impossible.

3. Give up on the "mpl_tookits" install path, and just install as "basemap". Advantages: easy, prevents further breakage when matplotlib changes the mpl_toolkits hierarchy. Disadvantage: Breaks user code. Users will have to change

from mpl_toolkits.basemap import Basemap

to

from basemap import Basemap.

Right now I favor solution (1).

Basemap.warpimage gives a wrong result when using a Hammer projection. Compare the top panel (Mollweide projection) to the middle panel (Hammer projection) in attached image.

Test script available from https://gist.github.com/ycopin/5012663

In its current implementation, basemap.addcyclic actually does not work for 2D arrays and fixing the number of input arrays to 2 is also somewhat arbitrary.

How about changing it to something like this:

def addcyclic(a,axis=-1):
    """
    ``arrout, lonsout = addcyclic((arrin, lonsin),axis=-1)``
    adds cyclic (wraparound) points in longitude to one or several arrays, 
    (e.g. ``arrin`` and ``lonsin``),
    where ``axis`` sets the dimension longitude is in (default: right-most).
    """
    def _addcyclic(a):
        aT = np.swapaxes(a,0,axis)
        idx = np.append(np.arange(aT.shape[0]),0)
        return np.swapaxes(aT[idx],axis,0)
    if isinstance(a,list) or isinstance(a,tuple):
        return map(_addcyclic,a)
    else:
        return _addcyclic(a)

Best,
Jonas

PS: Please let me know if this is stuff for a pull request, then I will create one.

Pull request #138 changed the behaviour of basemap.addcyclic. The data and longitude arrays are now being treated the same, i.e. the first point being appended to the end, but previously the longitude array was handled by adding the difference between longitudes to the last longitude, which is in my mind the correct behaviour.

An example:

import numpy as np
from mpl_toolkits.basemap import addcyclic

lons = np.arange(0, 360, 60)
data = np.ones([5, len(lons)]) * np.arange(len(lons))
new_data, new_lons = addcyclic(data, lons)
print new_lons

Prior to #138 being merged this would result in [0, 60, 120, 180, 240, 300, 360] but after it results in [0, 60, 120, 180, 240, 300, 0].

Here is the link to the flow visualization that I like very much

http://www.glerl.noaa.gov/res/glcfs/currents/glcfs-currents-month.php?mon=02&year=2013

It would be cool to be able to do it on my data....

Basemap.warpimage does not work when in celestial coordinates (see bottom panel on attached image).

Test script available from https://gist.github.com/ycopin/5012663

The penultimate example at http://matplotlib.github.com/basemap/users/examples.html shows the great circle from London to New York with a continent fill. Unfortunately (for those from mainland Europe) not all of the continents have been drawn.

My guess is that this is some data clipping issue related to the bounding box uses for clipping the GSHHS dataset.

Hi:

I was looking at the examples page here and noticed this artifact in the plot. See the attached image. Not sure if it is fixed in the latest version...

Sasha

I can't find a single direct link to the tarballs within the first couple of pages. This should be fairly prominent. I suggest having one in the sidebar, in the ToC on the front page (above Installing), and then also a link to it where it is mentioned to download the source tarballs in the install instructions.

Hey, It just occurred to me that my question I posted on StackOverflow might actually be a bug. So I'm going to put a ticket in here and hope that maybe someone can look and answer it.

http://stackoverflow.com/questions/13351157/create-simple-basemap-that-corresponds-to-a-specific-region-on-the-earth

The basic problem is that I can't create an image that is the proper size matching the projection/area desired.

If this is the wrong location for this, sorry in advance.

Thanks
Matt

fillcontinents() is not drawing all polygons that should be on the map.
Using the following code to visualize a map of Europe the main continent is not shown:

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
bm = Basemap(llcrnrlon=-10, llcrnrlat=35, urcrnrlon=25, urcrnrlat=60, projection='merc', resolution="c")
bm.fillcontinents()
plt.show()

Using a wider box in the basemap and after that adjusting the axes limits manually, the missing polygons are displayed correctly:

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
bm = Basemap(llcrnrlon=-30, llcrnrlat=15, urcrnrlon=55, urcrnrlat=80, projection='merc', resolution="c")
bm.fillcontinents()
ax = plt.gca()
ax.set_xlim(bm(-10, 35)[0], bm(25, 60)[0])
ax.set_ylim(bm(-10, 35)[1], bm(25, 60)[1])
plt.show()

I am using matlotlib 1.1.0 and basemap 1.0.3 on python 2.6.6.

(I couldnt find any option to attach my figures here..)

I'm not sure when this started, but I noticed it today...

When plotting a gridded dataset with contour/contourf I'm getting what appears to be excess cropping on the top and bottom of the image. (There is some excess cropping on the sides, but it's considerably less.) See the examples below:

However, when I plot the same data on the same map background and use pcolormesh, I don't have the same issues.

These images were created using the m.contourf(...), m.contour(...), and m.pcolormesh(...) methods. When I use straight matplotlib [plt.contourf(...) and plt.contour(...)], the images are correct (below)

I'm using master and have confirmed this behavior on both OSX and RHEL5.

In the docs of drawparallels it says that any other **kwargs are supplied to ax.text. The problem is that ax.text has a color parameter which I need to set somehow. But color is also a regular keyword of drawparallels which means that it won't get there.

I really, really, really don't want to use something like mpl.rcParams['text.color'] = colors.textcolor.

matplotlib-basemap-68a32f1/geos-3.3.3 $ make
[...]
mv -f .deps/Assert.Tpo .deps/Assert.Plo
mv -f .deps/Profiler.Tpo .deps/Profiler.Plo
mv -f .deps/GeometricShapeFactory.Tpo .deps/GeometricShapeFactory.Plo
/bin/sh ../../libtool --tag=CXX --mode=link g++ -DGEOS_INLINE -pedantic -Wall -ansi -Wno-long-long -ffloat-store -g -O2 -o libutil.la Assert.lo GeometricShapeFactory.lo math.lo Profiler.lo
libtool: link: ar cru .libs/libutil.a .libs/Assert.o .libs/GeometricShapeFactory.o .libs/math.o .libs/Profiler.o
libtool: link: ranlib .libs/libutil.a
libtool: link: ( cd ".libs" && rm -f "libutil.la" && ln -s "../libutil.la" "libutil.la" )
/bin/sh ../libtool --tag=CXX --mode=compile g++ -DHAVE_CONFIG_H -I. -I../include -I../include/geos -I../include -DGEOS_INLINE -pedantic -Wall -ansi -Wno-long-long -ffloat-store -g -O2 -MT inlines.lo -MD -MP -MF .deps/inlines.Tpo -c -o inlines.lo inlines.cpp
libtool: compile: g++ -DHAVE_CONFIG_H -I. -I../include -I../include/geos -I../include -DGEOS_INLINE -pedantic -Wall -ansi -Wno-long-long -ffloat-store -g -O2 -MT inlines.lo -MD -MP -MF .deps/inlines.Tpo -c inlines.cpp -fno-common -DPIC -o .libs/inlines.o
libtool: compile: g++ -DHAVE_CONFIG_H -I. -I../include -I../include/geos -I../include -DGEOS_INLINE -pedantic -Wall -ansi -Wno-long-long -ffloat-store -g -O2 -MT inlines.lo -MD -MP -MF .deps/inlines.Tpo -c inlines.cpp -o inlines.o >/dev/null 2>&1
mv -f .deps/inlines.Tpo .deps/inlines.Plo
/bin/sh ../libtool --tag=CXX --mode=link g++ -DGEOS_INLINE -pedantic -Wall -ansi -Wno-long-long -ffloat-store -g -O2 -release 3.3.3 -no-undefined -o libgeos.la -rpath /tools/geos/lib inlines.lo algorithm/libalgorithm.la geom/libgeom.la geomgraph/libgeomgraph.la index/libindex.la io/libio.la linearref/liblinearref.la noding/libnoding.la operation/liboperation.la planargraph/libplanargraph.la precision/libprecision.la simplify/libsimplify.la util/libutil.la
libtool: link: g++ -dynamiclib -o .libs/libgeos-3.3.3.dylib .libs/inlines.o -Wl,-force_load,algorithm/.libs/libalgorithm.a -Wl,-force_load,geom/.libs/libgeom.a -Wl,-force_load,geomgraph/.libs/libgeomgraph.a -Wl,-force_load,index/.libs/libindex.a -Wl,-force_load,io/.libs/libio.a -Wl,-force_load,linearref/.libs/liblinearref.a -Wl,-force_load,noding/.libs/libnoding.a -Wl,-force_load,operation/.libs/liboperation.a -Wl,-force_load,planargraph/.libs/libplanargraph.a -Wl,-force_load,precision/.libs/libprecision.a -Wl,-force_load,simplify/.libs/libsimplify.a -Wl,-force_load,util/.libs/libutil.a -O2 -install_name /tools/geos/lib/libgeos-3.3.3.dylib -Wl,-single_module
Undefined symbols for architecture x86_64:
"__ZNSt8auto_ptrIN4geos4geom8EnvelopeEEcvSt12auto_ptr_refIT_EIS2_EEv", referenced from:
virtual thunk to geos::geom::GeometryCollection::computeEnvelopeInternal() constin libgeom.a(GeometryCollection.o)
"std::auto_ptrgeos::geom::Envelope::auto_ptr(std::auto_ptr_refgeos::geom::Envelope)", referenced from:
virtual thunk to geos::geom::GeometryCollection::computeEnvelopeInternal() constin libgeom.a(GeometryCollection.o)
ld: symbol(s) not found for architecture x86_64
collect2: ld returned 1 exit status
make[2]: *** [libgeos.la] Error 1
make[1]: *** [all-recursive] Error 1
make: *** [all-recursive] Error 1

any idea?

TIA
dave

Kind of ran into an unexpected quirk with pcolormesh that perhaps could be improved to save users like me from themselves!

The X, Y arrays to this function need to be 1 larger than the data array so that the bounds of the grid cell are known. Matplotlib's code silently trims the array and will drop the rightmost column and top row.

The examples/plotsst.py suffers from this issue as X,Y,data are all (720, 1440), so the plot is not complete.

I am not sure what to suggest for a fix to this, perhaps simply documenting it is sufficient or maybe basemap should complain when this happens? If you think I should take my whining upstream to matplotlib, that's understandable :)

thank you.

Hi,

(Tested with basemap 1.0.6)

This works:

m.readshapefile('foo.shp', 'my_shape')

while this fails:

m.readshapefile('bar_0.1.shp', 'my_shape')

raise IOError('error reading shapefile %s.shp' % shapefile)
  fields = shf.fields
  coords = []; attributes = []

IOError: error reading shapefile bar_0.1.shp

Note that foo.shp and bar_0.1.shp are exactly the same file except for the filename. So it seems that the additional . causes the problem. Minimal example to reproduce behavior:

Rename huralll020.{shp,dfb,shx} to huralll020_0.1.{shp,dbf,shx} and change filename in examples/hurrtracks.py accordingly. Run hurrtracks.py.

Of course the obvious workaround is not to use . in filenames. But this is somewhat impracticable as I have a large number or files and would like keeping the filenames consistent accross filetypes (netCDF, shapfiles, pngs, etc).

Is there an easy fix?

Thanks very much,

Andy

Here is the example of the issue I am talking about (shown in the figure), the parallel intersects itself. The code that reproduces the image is here:
https://github.com/guziy/RPN/blob/master/src/demo/demo_basemap_problem.py

It reads the netcdf file from here:
https://github.com/guziy/RPN/blob/master/grid_arctic.nc

Thanks in advance for looking at this or for some insights on how to correct it.

Oleksandr (Sasha) Huziy

The script below generates a straight line in addition to the expected curved parallel. Although I have a patch that fixes this, it causes an error in the geos demos, so it is not a general solution.

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

m = Basemap(llcrnrlon=-178.5, llcrnrlat=84.5,
            urcrnrlon=-151.8, urcrnrlat=86.1,
            lon_0=-165, lat_0=85,
            projection='lcc')
m.drawparallels([85.2])
plt.show()

To make the environmental variable BASEMAPDATA function also for pyproj, please update pyproj.py with something like the following:

if 'BASEMAPDATA' in os.environ:
    pyproj_datadir = os.environ['BASEMAPDATA']
else:
    pyproj_datadir = os.sep.join([os.path.dirname(__file__), 'data'])

At the moment only the last line is present which causes problems when using py2exe

Hi, I asked out this on the email list and think I have resolved it

http://matplotlib.1069221.n5.nabble.com/RPM-conflicts-between-basemap-and-matplotlib-mpl-toolkits-init-py-oc-td39635.html

I think basemap should not own lib/mpl_toolkits/init.py in its MANIFEST.in as matplotlib already owns this file in its MANIFEST.in . I may be wrong, thanks for your time and wonderful library.

Jeff, I'm finishing wmsimage() function to Basemap class similar to arcgisimage() function and now testing.

Most WMS servers accept just epsg:4326 reference, but I can't figure how to reference coordinates. I would have expected that entering lat/lot pairs should do it, but I guess not. For example:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

plt.figure(figsize=(6,8))
epsg = 4326; lon1 = -85; lat1 = 24; lon2 = -79; lat2 =32
m=Basemap(epsg=epsg,resolution='i',llcrnrlon=lon1,llcrnrlat=lat1,\
           urcrnrlon=lon2,urcrnrlat=lat2)
m.arcgisimage(service='NatGeo_World_Map',verbose=True,xpixels=600)

returns this: bbox=-1,0,-1,0

or here is verbose url: http://server.arcgisonline.com/ArcGIS/rest/services/NatGeo_World_Map/MapServer/export?bbox=-1,0,-1,0&bboxSR=4326&imageSR=4326&size=600,800&dpi=96&format=png32&f=image

Why?
Or how can I use this projection in defined lat/lon pairs?

Hi,

I ran a code that is an example of "Plotting Maps: Visualizing Haiti Earthquake Crisis Data" on a book, Python for Data Analysis. Page 242-246

The code is supposed to create a plot map but I got an error with basemap module:

Traceback (most recent call last):
  File "Haiti.py", line 74, in <module>
    x, y = m(cat_data.LONGITUDE, cat_data.LATITUDE)
  File "/usr/local/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py", line 1148, in __call__
    xout,yout = self.projtran(x,y,inverse=inverse)
  File "/usr/local/lib/python2.7/site-packages/mpl_toolkits/basemap/proj.py", line 286, in __call__
    outx,outy = self._proj4(x, y, inverse=inverse)
  File "/usr/local/lib/python2.7/site-packages/mpl_toolkits/basemap/pyproj.py", line 388, in __call__
    _proj.Proj._fwd(self, inx, iny, radians=radians, errcheck=errcheck)
  File "_proj.pyx", line 122, in _proj.Proj._fwd (src/_proj.c:1571)
RuntimeError

I checked if basemap and proj module were installed okay on my Mac. Basemap was installed from source as instructed and proj was installed by Homebrew and they looks fine to me.

Any advice please?

Download code below

Haiti.csv file can be downloaded from https://github.com/pydata/pydata-book/raw/master/ch08/Haiti.csv

import pandas as pd
import numpy as np
from pandas import DataFrame

data = pd.read_csv('Haiti.csv')

data = data[(data.LATITUDE > 18) & (data.LATITUDE < 20) &
        (data.LONGITUDE > -75) & (data.LONGITUDE < -70)
        & data.CATEGORY.notnull()]

def to_cat_list(catstr):
    stripped = (x.strip() for x in catstr.split(','))
    return [x for x in stripped if x]

def get_all_categories(cat_series):
    cat_sets = (set(to_cat_list(x)) for x in cat_series) 
    return sorted(set.union(*cat_sets))

def get_english(cat):
    code, names = cat.split('.') 
    if '|' in names:
        names = names.split(' | ')[1] 
    return code, names.strip()

all_cats = get_all_categories(data.CATEGORY)
english_mapping = dict(get_english(x) for x in all_cats)

def get_code(seq):
    return [x.split('.')[0] for x in seq if x]

all_codes = get_code(all_cats)
code_index = pd.Index(np.unique(all_codes))
dummy_frame = DataFrame(np.zeros((len(data), len(code_index))),
                        index=data.index, columns=code_index)

for row, cat in zip(data.index, data.CATEGORY): 
    codes = get_code(to_cat_list(cat)) 
    dummy_frame.ix[row, codes] = 1

data = data.join(dummy_frame.add_prefix('category_'))

from mpl_toolkits.basemap import Basemap 
import matplotlib.pyplot as plt

def basic_haiti_map(ax=None, lllat=17.25, urlat=20.25, lllon=-75, urlon=-71):
    # create polar stereographic Basemap instance. 
    m = Basemap(ax=ax, projection='stere', 
                lon_0=(urlon + lllon) / 2, 
                lat_0=(urlat + lllat) / 2,
                llcrnrlat=lllat, urcrnrlat=urlat, 
                llcrnrlon=lllon, urcrnrlon=urlon, 
                resolution='f')
    # draw coastlines, state and country boundaries, edge of map. m.drawcoastlines()
    m.drawstates()
    m.drawcountries()
    return m

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 10)) 
fig.subplots_adjust(hspace=0.05, wspace=0.05)

to_plot = ['2a', '1', '3c', '7a']

lllat=17.25; urlat=20.25; lllon=-75; urlon=-71

for code, ax in zip(to_plot, axes.flat):
    m = basic_haiti_map(ax, lllat=lllat, urlat=urlat,
                        lllon=lllon, urlon=urlon) 

    cat_data = data[data['category_%s' % code] == 1]

    # compute map proj coordinates.
    print cat_data.LONGITUDE, cat_data.LATITUDE
    x, y = m(cat_data.LONGITUDE, cat_data.LATITUDE)

    m.plot(x, y, 'k.', alpha=0.5)
    ax.set_title('%s: %s' % (code, english_mapping[code]))

Now that I played with 'rotpole' a bit... I also found that dralsmask fails with ZeroDivisionError. This is the full traceback:

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-79-0ba416c748fa> in <module>()
----> 1 m.drawlsmask(resolution='c')

/home/razvan/.local/lib/python3.4/site-packages/mpl_toolkits/basemap/__init__.py in drawlsmask(self, land_color, ocean_color, lsmask, lsmask_lons, lsmask_lats, lakes, resolution, grid, **kwargs)
   3934         # to a rectangular map projection grid.
   3935             mask,x,y = self.transform_scalar(lsmask,lsmask_lons,-> 3936                        lsmask_lats,nx,ny,returnxy=True,order=0,masked=255)
   3937             lsmask_lats.dtype
   3938             # for these projections, points outside the projection

/home/razvan/.local/lib/python3.4/site-packages/mpl_toolkits/basemap/__init__.py in transform_scalar(self, datin, lons, lats, nx, ny, returnxy, checkbounds, order, masked)
   2961                 raise ValueError('grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)')
   2962         if returnxy:
-> 2963             lonsout, latsout, x, y = self.makegrid(nx,ny,returnxy=True)
   2964         else:
   2965             lonsout, latsout = self.makegrid(nx,ny)

/home/razvan/.local/lib/python3.4/site-packages/mpl_toolkits/basemap/__init__.py in makegrid(self, nx, ny, returnxy)
   1172         If ``returnxy = True``, the x,y values of the grid are returned also.
   1173         """
-> 1174         return self.projtran.makegrid(nx,ny,returnxy=returnxy)
   1175 
   1176     def _readboundarydata(self,name,as_polygons=False):

/home/razvan/.local/lib/python3.4/site-packages/mpl_toolkits/basemap/proj.py in makegrid(self, nx, ny, returnxy)
    332         if returnxy=True, the x,y values of the grid are returned also.
    333         """
--> 334         dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
    335         dy = (self.urcrnry-self.llcrnry)/(ny-1)
    336         x = self.llcrnrx+dx*np.indices((ny,nx),np.float32)[1,:,:]

ZeroDivisionError: float division by zero

I figured it's because of this:

nx = int((self.xmax-self.xmin)/dx)+1; ny = int((self.ymax-self.ymin)/dx)+1

here self.xmax, self.xmin, ... are used incorrectly for 'rotpole' giving nx == ny == 1.

I tried changing line 3928:

if self.projection == 'cyl':

to

if self.projection == 'cyl' or self.projection == 'rotpole':

with the effect of producing a nice pattern :) but no land sea mask unfortunately... so the issue is a bit dipper, but I don't have time to investigate now...

Edit: on version 1.0.7

I exported the GEOS_DIR like below

export GEOS_DIR = /Users/SurabhiAkshara/projects/basemap/geos-3.3.3/capi

and when i run python setup.py install , I get the following error. I even tried uninstalling Xcode, still facing the same error. Pls advise.

clang: error: unknown argument: '-mno-fused-madd' [-Wunused-command-line-argument-hard-error-in-future]
clang: note: this will be a hard error (cannot be downgraded to a warning) in the future
clang: error: unknown argument: '-mno-fused-madd' [-Wunused-command-line-argument-hard-error-in-future]
clang: note: this will be a hard error (cannot be downgraded to a warning) in the future
Traceback (most recent call last):
File "setup.py", line 104, in
objects = cc.compile(['nad2bin.c', 'src/pj_malloc.c'])
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/distutils/ccompiler.py", line 203, in CCompiler_compile
self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/distutils/unixccompiler.py", line 38, in UnixCCompiler__compile
raise CompileError(msg)
distutils.errors.CompileError: Command "cc -fno-strict-aliasing -fno-common -dynamic -arch x86_64 -arch i386 -g -Os -pipe -fno-common -fno-strict-aliasing -fwrapv -mno-fused-madd -DENABLE_DTRACE -DMACOSX -DNDEBUG -Wall -Wstrict-prototypes -Wshorten-64-to-32 -DNDEBUG -g -fwrapv -Os -Wall -Wstrict-prototypes -DENABLE_DTRACE -arch x86_64 -arch i386 -pipe -Isrc -c nad2bin.c -o nad2bin.o" failed with exit status 1

I have created two polar stereographic basemaps that look like this, sharing a plt.subplots figure:

http://imgur.com/ayqINQ6

When zooming in, this happens:

http://imgur.com/ca5OFgj

I have confirmed that this happens only for the Basemap with round=True option.

MPL: 1.2.1
basemap: 1.0.6
Backend: Qt4Agg

>>> Basemap(width=1000000, height=1000000, projection='stere', lat_0=-88, lon_0=-3)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "..\site-packages\mpl_toolkits\basemap\__init__.py", line 1085, in __init__
    self._readboundarydata('gshhs',as_polygons=True)
  File "..\site-packages\mpl_toolkits\basemap\__init__.py", line 1428, in _readboundarydata
    if not poly.is_valid(): poly=poly.fix()
  File "_geoslib.pyx", line 234, in _geoslib.BaseGeometry.fix (src/_geoslib.c:2233)
  File "_geoslib.pyx", line 334, in _geoslib.Polygon.__init__ (src/_geoslib.c:3479)
IndexError: index -1 is out of bounds for axis 0 with size 0

Strange thing is that e.g. using lon_0=-2 works.

I'm using basemap 1.0.8.

Any chance you'd consider supporting the Kavrayskiy VII, Eckert IV or Winkel-Tripel projections? It looks like PROJ.4 supports those projections.

Thanks!

There are small documentation gaps in Basemap 1.0.7. Specifically, the parameters o_lat_p and o_lon_p are not documented for the rotpole projection. It would be great if this could be fixed!

Cheers!

Dear Jeff,

I've been using basemap/netCDF4python/griddata for a few years now, and must say that they are some of the best and most useful tools I've ever come across. Thank you very much for your great work!

There is a question, however, that's been puzzling me for a long time now.

My data usually covers Greenland (or parts of it), and is in a stereographic projection, either "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs " (EPSG 3413) or "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" (no EPSG code).

Now I would like to set up a basemap instance based on the above projection definition, but where the center of the plot; (lon_0, lat_0) in basemap terms, is given by the center of the actual data I have. And not by (lon_0, lat_0) as given in the projection definition. I don't know how to do this.

What works is: I set up a basemap instance with

m = Basemap(width=width,
            height=height,
            resolution='i',
            projection='stere',
            lat_0=lat_0,
            lon_0=lon_0)

where width, height, are taken from the data, and (lon_0, lat_0) is the center coordinate of the data. This basemap instance is a stereographic projection which is somewhat different to the data. Therefore I have to convert the data into the basemap projection

xx, yy = m(lons, lats)

where lons and lats are the longitudes and latitudes of my data points. Then I can plot it with

m.pcolormesh(xx, yy, data)

Fine. Except that the data doesn't properly fill the plot area anymore, because it had to be rotated.

So, in conclusion: Is there a way to set up a basemap instance with the exact projection information given above, but where the center of the plot corresponds to the center of the data (and not of the lon_0, lat_0 of the projection)?

The script and data to produce the attached figures can be downloaded from http://pism-docs.org/download/stere_plot.zip

Cheers,

Andy

Dr. Andy Aschwanden
Postdoctoral Fellow
Geophysical Institute
903 Koyukuk Drive
Fairbanks AK 99775
USA

email [email protected]
phone +1 (907) 474-7199
fax +1 (907) 474-7290

I was playing with the projections and it took me a while to understand what's going on... I need the 'rotpole' projection to convert some locations and plot some graphs. I have all the information: lon_0 = 10, o_lon_p = 10, o_lat_p = 10, but I don't have the real lat/lon, I have the rotated one so I'm forced to use llcrnrx, llcrnry....

I have a climate model that gives me the output in this rotated pole projection with rotated left lon = -19.5, rotated right lon = 21, rotated bottom lat = -25, rotated top lat = 20.5.

First take:

import mpl_toolkits.basemap as mb
kwargs = {
    'projection': 'rotpole',
    'llcrnrx': -19.5,
    'llcrnry': -25,
    'urcrnrx': 21,
    'urcrnry': 20.5,
    'o_lon_p': 10,
    'o_lat_p': 32.5,
    'lat_0': 10
}
m = mb.Basemap(**kwargs)

Now m.drawcoastlines() will give me a huge map that goes backwards on X axis... After a full day of investigation... I saw that m.llcrnrx' takes the value of360-19.5 = 340.5andm.urcrnrx = 21` so this clearly dosn't work...

Is this behaviour wanted? I mean, when building m = mb.Basemap(...) using llcrnrlon, llcrnrlat... you can safely pass negative values to llcrnrlon but when using llcrnrx you're forced to use (0, 360)...

Working version is:

import mpl_toolkits.basemap as mb
kwargs = {
    'projection': 'rotpole',
    'llcrnrx': 0,
    'llcrnry': -25,
    'urcrnrx': 21,
    'urcrnry': 40.5,
    'o_lon_p': 10,
    'o_lat_p': 32.5,
    'lat_0': -9.5
}
m = mb.Basemap(**kwargs)

So now I'm forced to do:

# x, y taken from climate model
xoffset = 19.5
lon, lat = m([x + xoffset], [y], inverse=True)

Even more confusing is the fact that you can use negative values in convert as in m(-1, 0, inverse=True) and you get the correct negative real longitude value...

Edit: using version 1.0.7; fixed 'urcrnry' misspell