Basemap Tutorial

Basemap is a great tool for creating maps using python in a simple way. It's a `matplotlib <http://matplotlib.org/>`_ extension, so it has got all its features to create data visualizations, and adds the geographical projections and some datasets to be able to plot coast lines, countries, and so on directly from the library.

Basemap has got `some documentation <http://matplotlib.org/basemap/index.html>`_, but some things are a bit more difficult to find. I started a readthedocs page to extend a little the original documentation and examples, but it grew a little, and now covers many of the basemap possibilities.

Some of the examples from the tutorial

The tutorial can be found at http://basemaptutorial.readthedocs.org/, and all the examples and its source code, at GitHub and it's available for sharing or being modified by adding the attribution.

The tutorial covers:

• An introductory chapter, showing the basic functionalities, and how to work with the included methods for reading shapefiles
• The explanation and an example for each of the basemap methods
• A section with more advanced examples, such as 3D maps, multiplots or WRF examples

 I would really appreciate some feedback, the comments are open!

Basemap raster clipping with a shapefile

 Basemap is a great library for mapping faster than other python options, but there are some usual things I couldn't find how to do. Clipping a raster using a shape is one of them. Here's how do I do it

The output

As usual, all the code can be found at GitHub

Getting some data

The example plots some elevation data, taken from the SRTM. After looking for some options, the easiest to work with was this one: http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

The shapefile will be the border of Andorra, taken from Natural Earth

The result is a little poor because the resolution is low, but works well for the example.

The script

from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy
import shapefile

fig = plt.figure()
ax = fig.add_subplot(111)

sf = shapefile.Reader("ne_10m_admin_0_countries")

for shape_rec in sf.shapeRecords():
    if shape_rec.record[3] == 'Andorra':
        vertices = []
        codes = []
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i+1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)


m = Basemap(llcrnrlon=1.4,
    llcrnrlat=42.4,
    urcrnrlon=1.77,
    urcrnrlat=42.7,
    resolution = None, 
    projection = 'cyl')

ds = gdal.Open('srtm_37_04.tif')
data = ds.ReadAsArray()

gt = ds.GetGeoTransform()
x = numpy.linspace(gt[0], gt[0] + gt[1] * data.shape[1], data.shape[1])
y = numpy.linspace(gt[3], gt[3] + gt[5] * data.shape[0], data.shape[0])

xx, yy = numpy.meshgrid(x, y)

cs = m.contourf(xx,yy,data,range(0, 3600, 200))

for contour in cs.collections:
        contour.set_clip_path(clip)

plt.show()

• I used the pyshp library for reading the shapefile, since Fiona and GDAL don't work well together, and OGR was longer
• Lines 14 to 27 create the path. A Matplotlib path is made by two arrays. One with the points (called vertices in the script), and the other with the functions for every point (called codes)
• In our case, only straight lines have to be used, so there will be a MOVETO to indicate the beginning of the polygon, many LINETO to create the segments and one CLOSEPOLY for closing it
• Of course, only the polygon for Andorra has to be used. I get it from the shapefile attributes
•  The prt array is for managing multipolygons, which is not the case, but the code will create correct clipping for multipolygons
• The path is created using the Path function, and then added to a PathPatch, to be able to use it as a closed polygon. Note the trasnform=ax.transData attribute. This assumes the polygon coordinates to be the ones used in the data (longitudes and latitudes in our case). More information here
• Next code lines draw the map as usual. I have used a latlon projection, so all the values for the raster and shapefile can be used directly. If the output raster was in an other projection, the shapefile coordinates should be appended to the path using the output projection (m(pts[j][0], pts[j][1]))
• The x and y coordinates are calculated from the GDAL geotransform, and then turned into a matrix using meshgrid
• The clipping itself is made in the lines 48 and 49. For each drawn element, the method set_clip_path is applied

Links

SRTM data download 
Matplotlib path tutorial
Clipping an image in Basemap
Understanding Basemap transformations

Clipping a raster with GDAL python

Shortest distance to a geometry in a specified direction using Python

Looking at this map, I wondered how to calculate which geometry in a set is the closest to a point in a given direction.  

Usually, the problem is finding the closest geometry in general, which is easy using the distance function, but I couldn't find a solution for this other.

So I put me this problem: Which is the closest country that I have at each direction, knowing my geographical coordinates?

All the source code is, as usual, at GitHub

  The algorithm

The main idea is:

1. Create an infinite line from the point towards the desired direction.
2. Calculate the difference geometry between the  line and each polygon
1. If the polygon and the line actually intersect, the result will be a multi-line. The first line length of the multi-line is the distance we are looking for

So this would be the initial situation:

 And the distance to the polygon 1 would be calculated as:

The main problem is how to calculate the difference between the two geometries, but fortunately, shapely comes with this function, so coding it is not so difficult:

from shapely.geometry import Polygon
from shapely.geometry import LineString
from math import cos
from math import sin
from math import pi

def closest_polygon(x, y, angle, polygons, dist = 10000):
  
    angle = angle * pi / 180.0
    line = LineString([(x, y), (x + dist * sin(angle), y + dist * cos(angle))])

    dist_min = None
    closest_polygon = None
    for i in range(len(polygons)):
        difference = line.difference(polygons[i])
        if difference.geom_type == 'MultiLineString':
            dist = list(difference.geoms)[0].length
            if dist_min is None or dist_min > dist:
                dist_min = dist
                closest_polygon = i
        
    
    
    return {'closest_polygon': closest_polygon, 'distance': dist_min}


if __name__ == '__main__':

    polygons = []
    polygons.append(Polygon([(4, 2), (4, 4), (6, 4), (6, 2)]))
    polygons.append(Polygon([(7, 2), (7, 4), (9, 4), (9, 2)]))
   
    
    print closest_polygon(3, 3, 90, polygons)

• The main section creates the two squares using shapely
• The closest_polygon function calculates the closest polygon and its distance:
• A LineString to the desired direction is calculated. The dist is in the units used by the polygons. An infinite line isn't possible, so the distance must be larger than the further
• For each of the polygons to analyze, the difference is calculated using the shapely difference method
• Then, if the line and the polygon intersect (and the line is long enough), a MultilineString will be the result of the difference operation. The first String in the MultilineString is the one that connects our point with the polygon. Its length is the distance we are looking for

The example schema, drawn with the script draw_closest.py

Calculating the closest country in each direction

After getting the formula for calculating the closest polygon, the next step would be using it for something. So:

Which country do I have in all directions?

To create the script, some things have to be considered:

1. The projection should be azimuthal equidistant so the distances can be compared in all the directions from the given point
2. I've used the BaseMap library to draw the maps. I find it a bit tricky to use, but the code will be shorter

The script is used this way:

usage: closest_country.py [-h] [-n num_angles] [-o out_file] [-wf zoom_factor]
                          lon lat

Creates a map with the closest country in each direction

positional arguments:
  lon              The point longitude
  lat              The point latitude

optional arguments:
  -h, --help       show this help message and exit
  -n num_angles    Number of angles
  -o out_file      Out file. If present, saves the file instead of showing it
  -wf zoom_factor  The width factor. Use it to zoom in and out. Use > 1 to
                   draw a bigger area, and <1 for a smaller one. By default is
                   1


For example:

python closest_country.py -n 100 -wf 2.0 5 41

The code has some functions, but the main one is draw_map:

def draw_map(self, num_angles = 360, width_factor = 1.0):

        #Create the map, with no countries
        self.map = Basemap(projection='aeqd',
                    lat_0=self.center_lat,lon_0=self.center_lon,resolution =None)
        #Iterate over all the angles:
        self.read_shape()
        results = {}
        distances = []
        for num in range(num_angles):
            angle = num * 360./num_angles
            closest, dist = self.closest_polygon(angle)
            if closest is not None:
                distances.append(dist)
                if (self.names[closest] in results) == False:
                    results[self.names[closest]] = []

                results[self.names[closest]].append(angle)
        
        #The map zoom is calculated here, 
        #taking the 90% of the distances to be drawn by default       
        width = width_factor * sorted(distances)[
                int(-1 * round(len(distances)/10.))]

        #Create the figure so a legend can be added
        plt.close()
        fig = plt.figure()
        ax = fig.add_subplot(111)
        cmap = plt.get_cmap('Paired')
        
       
        self.map = Basemap(projection='aeqd', width=width, height=width,
                    lat_0=self.center_lat,lon_0=self.center_lon,resolution =None)
        self.read_shape()
        
        #Fill background.
        self.map.drawmapboundary(fill_color='aqua')

        #Draw parallels and meridians to give some references
        self.map.drawparallels(range(-80, 100, 20))
        self.map.drawmeridians(range(-180, 200, 20))

           
        #Draw a black dot at the center.
        xpt, ypt = self.map(self.center_lon, self.center_lat)
        self.map.plot([xpt],[ypt],'ko')
    
        #Draw the sectors
        for i in range(len(results.keys())):
            for angle in results[results.keys()[i]]:
                anglerad = float(angle) * pi / 180.0
                anglerad2 = float(angle + 360./num_angles) * pi / 180.0
                polygon = Polygon([(xpt, ypt), (xpt + width * sin(anglerad), ypt + width * cos(anglerad)), (xpt + width * sin(anglerad2), ypt + width * cos(anglerad2))])
                patch2b = PolygonPatch(polygon, fc=cmap(float(i)/(len(results) - 1)), ec=cmap(float(i)/(len(results) - 1)), alpha=1., zorder=1)
                ax.add_patch(patch2b)
        

        #Draw the countries
        for polygon in self.polygons:
            patch2b = PolygonPatch(polygon, fc='#555555', ec='#787878', alpha=1., zorder=2)
            ax.add_patch(patch2b)

        #Draw the legend
        cmap = self.cmap_discretize(cmap, len(results.keys()))
        mappable = cm.ScalarMappable(cmap=cmap)
        mappable.set_array([])
        mappable.set_clim(0, len(results))
        colorbar = plt.colorbar(mappable, ticks= [x + 0.5 for x in range(len(results.keys()))])
        colorbar.ax.set_yticklabels(results.keys())

        plt.title('Closest country')

• The first steps are used to calculate  the closest country in each direction, storing the result in a dict. The distance is calculated using the closest_polygon method, explained in the previous section..
• The actual map size is then calculated, taking the distance where the 90% of the polygons will appear. The width_factor can change this, because some times the result is not pretty enough. Some times has to much zoom and some, too few. Note that the aeqd i.e., Azimuthal Equidistant projection is used, since is the one that makes the distances in all directions comparable.
• Next steps are to actually drawing the map
• The sectors (the colors indicating the closest country) are drawn using the Descartes library and it's PolygonPatch
• The legend needs to change the color map to a discrete color map. I used a function called cmap_discretize, found here, to do it
• The legend is created using the examples found in this cookbook

Some outputs:

Next steps: What's across the ocean

Well, my original idea was creating a map like this one, showing the closest country when you are at the beach. Given a point and a direction (east or west in the example), calculating the country is easy, and doing it for each point in the coast is easy too. The problem is that doing it automatic is far more difficult, since you have to know the best direction (not easy in many places like islands), which countries to take as the origin, etc.
An other good thing would be doing the same, but with d3js, since  the point position could become interactive. I found some libraries like shapely.js or  jsts, but I think that they still don't implement the difference operation that we need.

Links 

The post examples source code

If you’re on the beach, this map shows you what’s across the ocean: The map that made me think about this problem

A LinkedIn discussion that gave me some ideas

Basemap links

How to install Basemap (you can use a virtual environment to test it without installing it in the whole system). Be sure to have pip installed, and the python-dev package in case you are using Ubuntu. Some distributions have Basemap as a system package too.
How to create an Azimuthal equidistant map with Basemap - The Azimuthal Equidistant projection
Some simple and useful Basemap examples
Advanced Basemap tricks that helped me to add the legend and much more
How to discretize a color map
Descartes: Drawing polygons in Matplotlib