Fill countries in python basemap

Question

Hi I am trying to plot a map using pythons basemap with some countries filled in a certain colour.

Is there a quick and easy solution out there??

Solution 1

As has already been said by @unutbu, Thomas' post here is exactly what you are after.

Should you want to do this with Cartopy, the corresponding code (in v0.7) can be adapted from http://scitools.org.uk/cartopy/docs/latest/tutorials/using_the_shapereader.html slightly:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shpreader
import itertools
import numpy as np

shapename = 'admin_0_countries'
countries_shp = shpreader.natural_earth(resolution='110m',
                                        category='cultural', name=shapename)

# some nice "earthy" colors
earth_colors = np.array([(199, 233, 192),
                                (161, 217, 155),
                                (116, 196, 118),
                                (65, 171, 93),
                                (35, 139, 69),
                                ]) / 255.
earth_colors = itertools.cycle(earth_colors)



ax = plt.axes(projection=ccrs.PlateCarree())
for country in shpreader.Reader(countries_shp).records():
    print country.attributes['name_long'], earth_colors.next()
    ax.add_geometries(country.geometry, ccrs.PlateCarree(),
                      facecolor=earth_colors.next(),
                      label=country.attributes['name_long'])

plt.show()

Solution 2

Inspired by the answer from pelson, I post the solution I have. I will leave it up to you which works best, so I will not accept any answer at the moment.

#! /usr/bin/env python

import sys
import os
from pylab import *
from mpl_toolkits.basemap import Basemap
import matplotlib as mp

from shapelib import ShapeFile
import dbflib
from matplotlib.collections import LineCollection
from matplotlib import cm

def get_shapeData(shp,dbf):
  for npoly in range(shp.info()[0]):
    shpsegs = []
    shpinfo = []

    shp_object = shp.read_object(npoly)
    verts = shp_object.vertices()
    rings = len(verts)
    for ring in range(rings):
        if ring == 0:
            shapedict = dbf.read_record(npoly)
        name = shapedict["name_long"]
        continent = shapedict["continent"]
        lons, lats = zip(*verts[ring])
        if max(lons) > 721. or min(lons) < -721. or max(lats) > 91. or min(lats) < -91:
            raise ValueError,msg
        x, y = m(lons, lats)
        shpsegs.append(zip(x,y))
        shapedict['RINGNUM'] = ring+1
        shapedict['SHAPENUM'] = npoly+1
        shpinfo.append(shapedict)

    lines = LineCollection(shpsegs,antialiaseds=(1,))
    lines.set_facecolors(cm.jet(np.random.rand(1)))
    lines.set_edgecolors('k')
    lines.set_linewidth(0.3)
    ax.add_collection(lines)


if __name__=='__main__':

  f=figure(figsize=(10,10))
  ax = plt.subplot(111)
  m = Basemap(projection='merc',llcrnrlat=30,urcrnrlat=72,\
            llcrnrlon=-40,urcrnrlon=50,resolution='c')
  m.drawcountries(linewidth=0.1,color='w')

  sfile = 'ne_10m_admin_0_countries'

  shp = ShapeFile(sfile)
  dbf = dbflib.open(sfile)
  get_shapeData(shp,dbf)

  show()
  sys.exit(0)

This is the result

Here is my example how to fill in Albania in the correct colour (not very elegant I know ;)).

  #HACK for Albania
  shpsegs = []
  shpinfo = []

  shp_object = shp.read_object(9)
  verts = shp_object.vertices()
  rings = len(verts)
  for ring in range(rings):
      if ring == 0:
          shapedict = dbf.read_record(9)
      name = shapedict["name_long"]
      continent = shapedict["continent"]
      lons, lats = zip(*verts[ring])
      if max(lons) > 721. or min(lons) < -721. or max(lats) > 91. or min(lats) < -91:
          raise ValueError,msg
      x, y = m(lons, lats)
      shpsegs.append(zip(x,y))
      shapedict['RINGNUM'] = ring+1
      shapedict['SHAPENUM'] = npoly+1
      shpinfo.append(shapedict)
  lines = LineCollection(shpsegs,antialiaseds=(1,))
  if name == 'Albania':
    lines.set_facecolors('w')
  lines.set_edgecolors('k')
  lines.set_linewidth(0.3)
  ax.add_collection(lines)

It is important that you do this after you have done all the other shapes. Perhaps you can get rid of some part of this code, but as I said it was sufficient for me.

For my application I coloured contries by name or continent, therefore these lines:

    name = shapedict["name_long"]
    continent = shapedict["continent"]

The data used I got from this website: http://www.naturalearthdata.com/

Solution 3

Updating @pelson answer for Python 3:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shpreader
import itertools
import numpy as np

shapename = 'admin_0_countries'
countries_shp = shpreader.natural_earth(resolution='110m',
                                        category='cultural', name=shapename)

print(countries_shp)

# some nice "earthy" colors
earth_colors = np.array([(199, 233, 192),
                         (161, 217, 155),
                         (116, 196, 118),
                         (65, 171, 93),
                         (35, 139, 69),
                        ]) / 255
earth_colors = itertools.cycle(earth_colors)

ax = plt.axes(projection=ccrs.PlateCarree())

for country in shpreader.Reader(countries_shp).records():
    print(country.attributes['NAME_LONG'], next(earth_colors))
    ax.add_geometries(country.geometry, ccrs.PlateCarree(),
                      facecolor=next(earth_colors),
                      label=country.attributes['NAME_LONG'])

plt.show()

Comments

[unutbu]

Maybe useful: geophysique.be/2011/01/27/…

[Fran Borcic]

I beleive this helps: matplotlib.1069221.n5.nabble.com/…

[red_tiger]

Thanks for those comments, they where most helpful. I also found a site with free country data, which was just what I was looking for: http://www.naturalearthdata.com/

[pelson]

@red_tiger - you could answer your own question with a small code snippet and output?

[Taryn]

Please note that you should post the essential parts of the answer here, on this site, or your post risks being deleted See the FAQ where it mentions answers that are 'barely more than a link'. You may still include the link if you wish, but only as a 'reference'. The answer should stand on its own without needing the link.

[pelson]

Thanks @bluefeet - I can see why that would be the case. I've updated the answer to give some new information (without duplicating the original link, which I did not own the copyright on). Cheers,

[red_tiger]

Yes, actually the same happens to Armenia. I had to to a work arround, by explicitly filling these two countries in afterwards. The inquiery with the people from naturalearthdata was not conclusive and I did not follow this up once I fixed it for me

[Leo]

Calling shpreader.natural_earth gives me an http 404 not found error, it apparently tries to download it?

[pelson]

That's right. There are gigabytes of data available to you on Natural Earth - it wouldn't be practical to install them all. If you update your cartopy version to v0.11.2 the recent change in NaturalEarth's URL scheme will be reflected, and you will no longer get a 404. github.com/SciTools/cartopy/pull/472 relates.

[tommy.carstensen]

@red_tiger I have the same problem with Argentina and Angola. Can you post your solution to the "Albanian problem"? What did the folks at NaturalEarth say? Thanks.

[tommy.carstensen]

@red_tiger Does Albania appear if you change the alpha? For example lines.set_alpha(0.5)?

[red_tiger]

@tommy.carstensen see my updated answer for your first question. I opened a topic in the NaturalEarth data forum (naturalearthdata.com/forums/topic/armenia-under-water) but it was not conclusive. And I did not follow it up after I got this hack implemented. And for your last question. I have not tried this, but do let me know if that has a positive effect.

[dd.]

I got a TypeError: 'Polygon' object is not iterable error running in Python 3.8, which corresponded with (this GitHub issue)[github.com/SciTools/cartopy/issues/948]. I fixed it by making country.geometry into a list like this: ax.add_geometries([country.geometry], ccrs.PlateCarree(), ...

[Jan Bodnar]

@dd The error can be fixed with the following adjustement: ax.add_geometries([n.geometry], ccrs.PlateCarree()...