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" [email protected] 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 https://github.com/matplotlib/basemap/issues/165#issuecomment-69270271.

[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.

[jenshnielsen]

Another work around is to use the Transverse Mercator Projection http://matplotlib.org/basemap/users/tmerc.html

I the Belgium example that seems to return the expected result.

[Rabijns]

sorry, clicked a little to fast. It should indeed remain open....

[Rabijns]

My idea about lat0 and lon0 does not work.

I can confirm the info of jenshnielsen that with transverse Mercator projection, the scale is correct. The only change to the code is: projection='tmerc' instead of 'merc'

[Rabijns]

I played some more with the two different mercator projections. I also set the labelstyle to fancy, because the scale will then include more info. When the labelstyle is set to fancy, and when using the mercator projection, the scale does mention a 0.63 factor.

map.drawmapscale(4.27,51.11,4.5,50.89,228, barstyle='fancy', units='km', labelstyle='fancy', fontsize= '7')

[lafaysse]

For people who want to use that function with a trick to fix this bug: it is possible to modify the labels of the scale. Here is my code to plot a correct scale of length 50 km over French Alps:

   # Map definition
   latmin=43.9
   latmax=46.5
   lonmin=5.2
   lonmax=7.9        

   mymap = Basemap(llcrnrlat=latmin, urcrnrlat=latmax, llcrnrlon=lonmin, urcrnrlon=lonmax, projection="merc", resolution='i', fix_aspect=True)

    # Distance we want to plot
    dref=50

    # Coordinates
    lat0=mymap.llcrnrlat+0.2
    lon0=mymap.llcrnrlon+0.4

    # Tricked distance to provide to the the function 
    distance=dref/np.cos(lat0*np.pi/180.)

    # Due to the bug, the function will draw a bar of length dref
    scale=mymap.drawmapscale(lon0,lat0,lon0,lat0,distance,barstyle='fancy', units='km', labelstyle='simple',fillcolor1='w', fillcolor2='#555555', fontcolor='#555555')        

    # Modify the labels with dref instead of distance
    scale[12].set_text(dref/2)
    scale[13].set_text(dref)