24

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??

bmu
  • 35,119
  • 13
  • 91
  • 108
red_tiger
  • 1,402
  • 3
  • 16
  • 32
  • 2
    Maybe useful: http://www.geophysique.be/2011/01/27/matplotlib-basemap-tutorial-07-shapefiles-unleached/ – unutbu Nov 15 '12 at 12:03
  • I beleive this helps: http://matplotlib.1069221.n5.nabble.com/How-to-draw-a-specific-country-by-basemap-td15744.html – Fran Borcic Nov 15 '12 at 15:19
  • 2
    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/](http://www.naturalearthdata.com/) – red_tiger Nov 16 '12 at 09:27
  • 1
    @red_tiger - you could answer your own question with a small code snippet and output? – pelson Jan 04 '13 at 09:45

3 Answers3

20

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()

output

pelson
  • 21,252
  • 4
  • 92
  • 99
  • 5
    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'.](http://stackoverflow.com/faq#deletion) 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. – Taryn May 08 '13 at 21:12
  • 1
    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, – pelson May 08 '13 at 21:33
  • Calling shpreader.natural_earth gives me an http 404 not found error, it apparently tries to download it? – Leo Dec 05 '14 at 13:38
  • 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. https://github.com/SciTools/cartopy/pull/472 relates. – pelson Dec 05 '14 at 18:31
10

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

example for filling in countries in different colours

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/

red_tiger
  • 1,402
  • 3
  • 16
  • 32
  • 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 – red_tiger May 15 '13 at 13:44
  • @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 Jan 14 '16 at 17:36
  • @red_tiger Does Albania appear if you change the alpha? For example `lines.set_alpha(0.5)`? – tommy.carstensen Jan 15 '16 at 01:09
  • 1
    @tommy.carstensen see my updated answer for your first question. I opened a topic in the NaturalEarth data forum (http://www.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. – red_tiger Jan 15 '16 at 09:30
10

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()
user1718097
  • 4,090
  • 11
  • 48
  • 63
  • 1
    I got a ```TypeError: 'Polygon' object is not iterable``` error running in Python 3.8, which corresponded with (this GitHub issue)[https://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(), ...``` – dd. May 05 '20 at 10:57
  • 1
    @dd The error can be fixed with the following adjustement: `ax.add_geometries([n.geometry], ccrs.PlateCarree()...` – Jan Bodnar May 12 '20 at 16:07