Easiest way to plot data on country map with python

Question

Could not delete question. Please refer to question: Shade states of a country according to dictionary values with Basemap

I want to plot data (number of sick people for a certain year) on each state of Mexico. I am using jupyter notebook. So far I have seen several options and tutorials, but none seem to seem to explicitly explain how to plot the map of a country. Below I explain some options/tutorial I have seen and why they have not worked (this I do just to argue that tutorials are not very straight forward):

1. Bokeh (http://docs.bokeh.org/en/latest/docs/gallery/texas.html). In the tutorial texas state is plotted given that us_counties is in bokeh.sampledata. However I have not found other countries in the sampledata.

2. mpl_toolkits.basemap (http://www.geophysique.be/2011/01/27/matplotlib-basemap-tutorial-07-shapefiles-unleached/). Although I am able to import shapefile, I cannot run from shapefile import ShapeFile (ImportError: cannot import name ShapeFile). Furthermore I have not been able to download dbflib library.

3. Vincent (Why Python Vincent map visuzalization does not map data from Data Frame?) When I run the code from the answer in said tutorial no image appears (even though I used command vincent.core.initialize_notebook() ).

4. Plotly (https://plot.ly/python/choropleth-maps/). The tutorial plots the map of USA importing information from a csv table (no information of other countries available). If wanting to plot another country, would it be possible to make the table?

Explored this 4 options I have found tutorials not to be very clear or easy to follow. I find it hard to believe that plotting a map of a country is difficult in python. I think there must be an easier way than the ones explained in the past tutorials.

The question is: Which is the easiest (hopefully simple) way to plot the map of a certain country (any) with python and how?

I have installed the following packages: matplotlib, pyshp, mpl_toolkits.basemap, bokeh, pandas, numpy. I have also downloaded Mexico's map from http://www.gadm.org/

Thanks in advance.

Answer 1

this might not be the answer you hope for, but have you looked at Plotly for maps? Their example looks exactly like what you want to do, although I'm not familiar with it myself, and I don't know which maps they have available and how easy it would be to upload your own.

Answer 2

While this question seems to be unanswerable in its current form, I'll at least note that you seem to be something wrong when using basemap - you don't want to import Shapefile, but simply read it using the readshapefile method of a Basemap object like so:

m = Basemap(projection='tmerc')
m.readshapefile("/path/to/your/shapefile", "mexican_states")

You will then be able to access the coordinates of each state's boundaries via m.mexican_states (as a list of arrays) and the corresponding information (such as names, maybe an indentifying code) by m.mexican_states_info. You will then need some sort of dict or DataFrame containing names/codes for states (corresponding to what is in m.mexican_states_info) and the values you want to plot. A simple example would work something like this, assuming that you have a dict called mexican_states_sick_people that looks something like {"Mexico City":123, "Chiapas":35, ...}:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Polygon
from descartes import PolygonPatch

fig, ax = plt.subplots()

# Set up basemap and read in state shapefile (this will draw all state boundaries)
m = Basemap(projection='tmerc')
m.readshapefile("/path/to/your/shapefile", "mexican_states")

# Get maximum number of sick people to calculate shades for states based on relative number    
max_sick = np.max(mexican_states_sick_people.values())

# Loop through the states contained in shapefile, attaching a PolygonPatch for each of them with shade corresponding to relative number of sick people
state_patches = []
for coordinates, state in zip(m.mexican_states, m.mexican_states_info):
    if state["State_name"] in mexican_states_sick_people.keys():
        shade = mexican_states_sick_people[state["State_name"]]/max_sick       
        state_patches.append(PolygonPatch(Polygon(coordinates), fc = "darkred", ec='#555555', lw=.2, alpha=shade, zorder=4)

 # Put PatchCollection of states on the map
ax.add_collection(PatchCollection(state_patches, match_original=True))

This example should be more or less functional, if you have a working shapefile of states and make sure that the dataset of sick people you have has some sort of identifier (name or code) for each state that allows you to match up the numbers with the identifier of states in the shapefile (this is what the shade = ... line in the loop relies upon - in the example I'm accessing the vals in the dictionary using a name from the shapefile as key).

Hope this helps, good luck!