Soil cadmium pollution over Europe

make_map.py

{
 "metadata": {
  "name": "Cd2+ soil pollution in Europe"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "First the necessary imports and settings:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from os.path import expanduser, expandvars\n\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom mpl_toolkits.basemap import Basemap, cm",
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "The map data are taken from http://eusoils.jrc.ec.europa.eu/foregshmc/HM_coordinate_system.txt. The used reference system is EPSG3035, which was found out by putting the last few lines of the `HM_coordinate_system.txt` file (those between the last `>>` and `<<` tags) into [this website](http://prj2epsg.org/epsg/3035). The lats and lons of upper left/lower right corner were calculated using the following command:\n\n    > gdaltransform -s_srs epsg:3035 -t_srs \"+proj=latlong\"\n    2636000 1423000\n    -8.21935120696638 34.0515805814956 0\n    5956000 5418000\n    50.467681487278 67.3886909774455 0\n    \nThe first values are the longitudes, second ones are the latitudes, third values I don't fucking care.\nSo, first we create the map and load the data."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# This is the figure we want to draw the map to...\nf = plt.figure(edgecolor='white')\n# ...and this is the figure's plot\nax = f.add_subplot(111, frame_on=False)\n# Create the map with reference system EPSG-3035\nm = Basemap(epsg=3035,\n            # High resolution\n            resolution='l',\n            # lower left corner longitude/latitude\n            llcrnrlon=-8.21935120696638,\n            llcrnrlat=34.0515805814956,\n            # upper right corner longitude/latitude\n            urcrnrlon=50.467681487278,\n            urcrnrlat=67.3886909774455,\n            # Don't display map features smaller than 25 km^2\n            area_thresh=100,\n            # Draw the map into the \"ax\" plot\n            ax=ax,\n            )\n\n# Load the data from a file, ignore the first six lines within that file as they are headers\ndata = np.loadtxt(expanduser(expandvars('~/dls/foregsHMC5km/Cd5km.asc')), skiprows=6)\n# Create a grid of longitudes and latitudes on which we can plot the data\nlons, lats = m.makegrid(data.shape[1], data.shape[0])\n\n\n# Draw the continents in light gray\nm.fillcontinents(color='0.8', ax=ax)\n# Draw country borders in white\nm.drawcountries(linewidth=0.4, color='white', zorder=100, ax=ax)\n# Set color thresholds\nclevs = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]\n# Plot data onto map\ncs = m.contourf(lons, lats, data[::-1], clevs, latlon=True, zorder=99, cmap='Spectral_r', ax=ax)\n# Add legend\ncbar = m.colorbar(cs, location='left', pad=\"5%\", fig=f, ax=ax)\n# Remove legend outline\ncbar.ax.artists.remove(cbar.outline)\n# Add legend label\ncbar.set_label('mg Cd/kg soil')\n# Put legend label to the left\ncbar.ax.yaxis.set_label_position('left')\n# Put legend ticks to the left\ncbar.ax.yaxis.tick_left()",
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## References\n\n* http://eusoils.jrc.ec.europa.eu/foregshmc/\n* Rodriguez Lado, L., Hengl, T., Reuter, H.I., (2008) [Heavy metals in European soils: a geostatistical analysis of the FOREGS Geochemical database.](http://www.sciencedirect.com/science/article/pii/S0016706108002668) Geoderma 148, 189-199. "
    }
   ],
   "metadata": {}
  }
 ]
}