[GIS] Using Google Maps image as layer with Matplotlib Basemap

choroplethgoogle mapslayersmatplotlib-basemapshapefile

I have produced a heat map of assaults in Sydney by suburb. However, my map looks a bit plain and I'd like to do something fancier – I'd like to use a Google Maps image as a background layer – I haven't figured out how to use a Google Maps image as a layer just yet and I'd like your help. I would like something that looks like this: http://www.ncomva.se/guide/Images/Layers.PNG As you can see, a Google Maps image is used a base layer. Then on top of the Google Maps image they put a see through heat map. Basically, I want to know is if this is possible in Python – I've tried to figure out how to do this but whatever I've tried just isn't working. I'm new to Python so this is difficult! If this is possible, can you please direct me to code that would be helpful?

I've also included some code below (but not the full code as it's too long) that gives you an idea how to use basemap. 

import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import Normalize
from matplotlib.collections import PatchCollection
import matplotlib.artist
import matplotlib.colorbar
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Point, Polygon, MultiPoint, MultiPolygon, GeometryCollection
from pysal.esda.mapclassify import Natural_Breaks as nb
import pysal.esda.mapclassify as mapclassify
from descartes import PolygonPatch
import fiona
from itertools import chain
import numpy as np
from shapely.geometry import asShape, mapping

#### http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_excel.html
data = pd.ExcelFile('data/sr16-14330 Jane Baseby Monthly counts of non-DV assault.xls')
data = data.parse('Suburb', skiprows=4, skip_footer=5, index_col=None, na_values=['NA'])

# load in the shapefile - use the file SSC_2016_AUST or "suburb_subset" for only a subset of suburbs
shp = fiona.open('data/ssc_2016_aust_shape/suburb_subset.shp')
bds = shp.bounds  # extract shapefile boundary
shp.close()

extra = 0.1  # define a variable called extra which we will use for padding the map when we display it (in this case I've selected a 10% pad)
#ll = (bds[0], bds[1])  # define the lower left hand boundary (longitude, latitude)
#ur = (bds[2], bds[3])  # define the upper right hand boundary (longitude, latitude)

NSEW = {'S': -33.892078720829694, 'E': 151.24185030131247, 'W': 151.18691866068747, 'N': -33.846468362714106}

ll = (NSEW["W"],NSEW["S"])
ur = (NSEW["E"],NSEW["N"])

coords = list(chain(ll, ur))  # concatenate the lower left and upper right into a variable called coordinates
w, h = coords[2] - coords[0], coords[3] - coords[1]  # define variables for the width and the height of the map

m = Basemap(
    # set projection to 'tmerc' which is apparently less distorting when close-in
    projection='tmerc',

    # set longitude as average of lower, upper longitude bounds
    lon_0=np.average([bds[0], bds[2]]),

    # set latitude as average of lower,upper latitude bounds
    lat_0=np.average([bds[1], bds[3]]),

    ellps='WGS84',

    # set the map boundaries. Note that we use the extra variable to provide a 10% buffer around the map
    llcrnrlon=coords[0],
    llcrnrlat=coords[1],
    urcrnrlon=coords[2],
    urcrnrlat=coords[3],

    # resolution of boundary database to use. Can be c (crude), l (low), i (intermediate), h (high), f (full) or None.
    resolution='f',

    # alpha = 0.5,

    # don't show the axis ticks automatically
    suppress_ticks=True
    #,epsg=4269
)

m.readshapefile('data/ssc_2016_aust_shape/test', 'assaults', color='none', zorder=2)

Best Answer

I figured out how to put a satellite image as a background layer. All I needed to do was add in the line:

m.arcgisimage(service='ESRI_Imagery_World_2D', xpixels=2000, verbose=True, zorder=1

and make sure that basemap has "epsg=4326" as an argument (it basically lets basemap know to use Sydney as a background image)

Related Question