scale created by drawmapscale seems incorrect

[Rabijns]

Issue description

Drawing a scale in basemap using 'drawmapscale' results in a scale which is not correct. It seems the bar is about 1.6 times to short. So, a scale with the indication 10km above it is actually only about 6.25 km long.

The two attachments show:

• a basemap plot (code below) for a region around Belgium. Belgium is about 230 km width, when measured due east starting from the most southern point of the coastline. As can be seen, the ruler is to small, it does not go from coast till the eastern border.
• a screenshot from google earth showing the same area with a measurement, confirming that the country is indeed 230 km wide.

Code which shows the issue:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
map = Basemap(projection='merc', lat_0 = 50, lon_0 = 4,resolution = 'h', area_thresh = 0.1,
llcrnrlon=1.75, llcrnrlat=49.0,urcrnrlon=7, urcrnrlat=52)
map.drawcoastlines()
map.drawcountries()
map.drawmapboundary(fill_color = 'lightblue')
map.fillcontinents(color = 'white',lake_color = 'lightblue')
map.drawmapscale(4.27,51.11,4.5,50.89,228, barstyle='fancy', units='km', labelstyle='simple',\
fillcolor1='w', fillcolor2='#555555', fontcolor='#555555')
plt.show(block=False)

I tried similar plottings in nautical mile, meters, etc, all with the same problem: the ruler shows the correct numbers, but is 1.6x to small.
I also tried with much smaller scales (eg 10 km),other barstyle, etc and referenced them to known distances on a map, and the same issue occurs.

Technical info:

Python version: 2.7.6
matplotlib version: 1.3.1
Basemap 1.0.7 + dfsg-1
OS: linux Mint 17 Qiana (64-bit version)
Installation method: all packages installed automatically during Mint installation or via Synaptic package manager using the official repository
no customization

[Rabijns]

I did some further testing: the issue is also seen with basemap 1.0.5.

I also took a look at the basemap source code for drawmapscale.
It seems the basic idea for the scale is to go from lat/lon to x/y, add the required distance (variable "length") and go back to lat/lon. So, I recreated these steps (see code below), using two points with the same latitude. The conclusions seem to be:

• the distance determined by the tool between the two points (dist12) is in line with the distance reported by different websites
• going forward and back does not result in an identical longitude coordinate for the second point.
• in the drawmapscale code, a scalefactor is calculated (variable "scalefact"), but it is only displayed in a legend, not used elsewhere (the value is by the way not visible above the fancy legend in my screenshots above, so possibly this is added in a version later than 1.0.7).
• in my example below, the ratio between the distance detemined directly (dist12) and after going back and forward (dist12r) is 1.555729, which is very close to the estimated error of 1.6 in my original post.

Based on this, I wonder whether it is correct to simply add the requested length of the scale in x/y without any correction for the scale factor. I have however zero knowledge on projection and coordinate systems, so I simpy do not know.

One alternative solution could be to simply not go from lat/lon to x/y and back, but to do all calculations in lat/lon. The points of the legend could then be determined using pyproj functions fwd and npts.

from mpl_toolkits.basemap import Basemap, pyproj

# create two points and define projection
lon1 = 5. ; lat1 = 50.
lon2 = 6. ; lat2 = 50.

Proj = Basemap(projection='merc', lat_0 = lat1, lon_0 = lon1,
        resolution =  'h', area_thresh = 0.1,
        llcrnrlon=0.5*lon1, llcrnrlat=0.9*lat1,
        urcrnrlon=2.*lon2, urcrnrlat=1.1*lat2)

# distance calculation between points
gc = pyproj.Geod(a=Proj.rmajor,b=Proj.rminor)
az1, az2, dist12 = gc.inv(lon1,lat1,lon2,lat2)
print "distance based on Lat/lon coordiantes: %f" %(dist12)

# calculate x, y coordinates for the second point similar as in drawmapscale,
# by going to x/y, adding distance and go back to lat/lon
x1,y1 = Proj(lon1,lat1)
lon2_reverse, lat2_reverse = Proj(x1+dist12,y1,inverse = True)
print "original Lat/lon for second point: %f, %f \n calculated lat/lon: %f, %f" \
        %(lat2, lon2, lat2_reverse, lon2_reverse)

# determine scale factor
az1, az2, dist12r = gc.inv(lon1,lat1,lon2_reverse,lat2_reverse)
print "ratio of distance: %f" %(dist12/dist12r)

Resulting output:
distance based on Lat/lon coordiantes: 71474.155087
original Lat/lon for second point: 50.000000, 6.000000
calculated lat/lon: 50.000000, 5.642783
ratio of distance: 1.555729

[mankoff]

This seems to be a fairly serious bug. I have a different interpretation for its existence. I think it is due to the mercator projection. The distance is correct when plotted at (0,0), where 1 degree is 111 km (longitude or latitude, roughly). If you move north, it maintains 111 km per degree longitude, even though that is not correct, since longitude lines approach each other.

Regardless, I saw another issue on this site mentioning that Basemap is being deprecated in favor of Cartopy? So this might not get fixed...

[mankoff]

Cartopy doesn't have a scalebar function, but here is a 3rd party function that will create one. http://nbviewer.ipython.org/github/pp-mo/iris_example_code/blob/cartopy_scalebar/map_scalebar.ipynb

[WeatherGod]

Basemap isn't going away any time soon, but will likely not get new
features. Serious bugs like this should get fixed, and your observation
seems to be an important clue.
On Jan 8, 2015 6:36 PM, "Ken Mankoff" notifications@github.com wrote:

Cartopy doesn't have a scalebar function, but here is a 3rd party function
that will create one.

—
Reply to this email directly or view it on GitHub
#165 (comment).

[efiring]

Chiming in: Basemap is a wonderful and crucial package--many thanks to @jswhit for the enormous amount of work he has put into it over the years. I expect it will be well worth maintaining for quite some time. Cartopy looks very nice, too, but it appears to be quite a long way away from having all the characteristics that would enable me to use it as a replacement for Basemap.

[mankoff]

Here is a code snippet demonstrating the bug which appears in merc but not in ortho. Latitude degree is 111 km at the equator. I'm making the scale 30*2 degrees wide. The northern bar in ortho appears wrong, but I don't think it is. The scale is based at the center of the bar, and it doesn't take into account the curvature of the earth elsewhere where the bar is drawn, nor should it.

from mpl_toolkits.basemap import Basemap

fig = plt.figure('scale', figsize=(4,3))
fig.clf()
fig.set_tight_layout(True)

ax = fig.add_subplot(111)

lat, lon = 0, 0
m = Basemap(projection='merc',
            lon_0 = lon, lat_0 = lat,
            llcrnrlon=lon-45,
            urcrnrlon=lon+45,
            llcrnrlat=-80,
            urcrnrlat=80,
            resolution='c')

# m = Basemap(projection='ortho',
#             lon_0 = lon, lat_0 = lat,
#             resolution='c')

m.drawparallels(np.arange(-80,80,10), color='gray')  # 10 deg lat bands
m.drawmeridians(np.arange(-90,90,30), color='gray')  # 30 deg lon bands
m.drawcoastlines(color='gray')

m.drawmapscale(0,0, 0,0, 111*2*30, barstyle='fancy', units='km')
m.drawmapscale(0,45,0,45, 111*2*30, barstyle='fancy', units='km')

[mankoff]

Perhaps related, the tissot function exhibits a similar bug.

Add the following two lines to the above code:

m.tissot(0,0,30,90,alpha=0.5)
m.tissot(0,45,30,90,alpha=0.5)

This should draw a circle with a radius of 30 degrees at the equator, and again at 45 degrees north.

The two screenshots show the 45 degree north circle is not 30 degrees (note: For the ortho projection, I've rotated it so 45 degrees is the center latitude by setting lat_0=45).

[mankoff]

Last comment for a while, in reply to @Rabijns (OP). You say scale is off by 1.6. Your center latitude is 50. Note that np.cos(50*(pi/180.)) = 0.64, and 1/0.64 = 1.55.

[Rabijns]

Thanks to all of you for looking into my problem.
@ Mankoff: this 1.55 seems to be close or identical to the 1.555 I mentioned in my second post.

I wonder if a workaround could be to tell drawmapscale to do its calculation at the equator (via the lat0 and Lon0 parameters)
From the documentation: Draw a map scale at lon,lat of length length representing distance in the map projection coordinates at lon0,lat0.
So tyhere is some kind of notion that the scale differs depending on where you calculate its length.

I'll try and post the results.

As a temporary work-around a wrote my own function to draw a mapscale. For this, I calculate the postion of the scale bars directly in lat/lon via range/bearing and not via the projection coordinate system.

[mankoff]

This bug report should remain open.