By Shaheen Syed and Charlotte Weber,
This example shows how to plot and color European seas with Python, Basemap and the ICES shape file. It serves as a basic example from which more elements such as annotations and other plotting elements can be added.
It might be the case that during your research you found some fancy numbers that relate to the different seas within Europe. Let’s say you have data on stocks, landings, vessels or some other valuable number that you would like to visually represent. The basemap toolkit is one way of visually representing your data onto a map. It’s part of Matplotlib, the de facto standard when creating plots with Python. To use basemap we need something that indicates the countries, regions, states etc so we can use them within our plots. That’s where shapefiles come in. Without going into too much detail, the shapefile makes sure we can pinpoint these countries, regions, states within our map. But what about the seas? We’ll have to thank ICES here.
Shapefile for ecoregions (European seas) are provided by ICES and can be viewed and downloaded here.
How to use them? We’ve created a simple example with some dummy data to color the regions based on some arbitrary frequency. Of course, you can easily load data from a CSV file and have the map say something meaningful. Hope this helps.
# imports
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib import cm
import matplotlib
import numpy as np
import shapefile
##################################
####### FUNCTIONS ########
##################################
def get_color(count):
cmap = matplotlib.cm.get_cmap('Blues')
# we normalize the scale here, adjust values accordingly
norm = matplotlib.colors.Normalize(vmin=0.0, vmax=300.0)
return cmap(norm(count))
##################################
####### MAIN #######
##################################
# define figure and axes
fig, axs = plt.subplots(1, 1, figsize=(10,10))
# define coordinates for map
# llcrnrlon longitude of lower left hand corner of the desired map domain (degrees).
# llcrnrlat latitude of lower left hand corner of the desired map domain (degrees).
# urcrnrlon longitude of upper right hand corner of the desired map domain (degrees).
# urcrnrlat latitude of upper right hand corner of the desired map domain (degrees).
x1 = -30.
x2 = 50.
y1 = 32.
y2 = 70.
# create map
m = Basemap(resolution='i', projection='merc', llcrnrlat=y1, urcrnrlat=y2, llcrnrlon=x1, urcrnrlon=x2)
# draw bounderies on map
m.drawcountries(linewidth=0.5)
m.drawcoastlines(linewidth=0.5)
# data of seas [0]=label (not used within map), [1]=frequency of some value, [2]=id for sea within shape file
# this list of list is part of the code but can easily be loaded from e.g a csv file
# note that we are not using all shape ids
data_seas_rows = [
['Greenland Sea', 28, 1 ],
['Bay of Biscay & Iberian Peninsula', 90, 2 ],
['Azores', 200, 3 ],
['western Mediterranean', 200, 4 ],
['Ionian Sea', 104, 5 ],
['Black Sea', 100, 6 ],
['Adriatic', 120, 7 ],
['Aegean Sea', 300, 8 ],
['Celtic Sea', 10, 9 ],
['Baltic Sea', 170, 10 ],
['North Sea', 58, 11 ],
['Iceland Sea', 20, 13 ],
['Faroe Islands', 23, 15 ],
['Norwegian Sea', 190, 16 ],
['North East Atlantic', 79, 17 ]
]
# Load shapefile from file
r = shapefile.Reader(r"ICES-Shape/ICES_ecoregions_20150113_no_land")
# plot every row from data_seas_rows
for row in data_seas_rows:
# get the shape id from our list, which is listed in position 2 (note the zero-based counting)
i = row[2]
# now point to the right shape and record from the ICES shape file
# note that we start with position i-1 to i
shapes = r.shapes()[i-1:i]
records = r.records()[i-1:i]
# use records and shapes
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_edgecolors('black')
lines.set_linewidth(0.1)
# color the shape of the sea
# we call function get_color here and have as argument the 2nd item from our data_seas_rows list
# basically this can be anything that you study e.g. frequency of catch, vessels etc.
lines.set_facecolors(get_color(row[1]))
# add to axes
axs.add_collection(lines)
# save the plot
plt.tight_layout()
plt.savefig('map.png', dpi=300)