This example is taken out of my Ph.D. thesis. It shows the North Sea bathymetry and the topography of the surrounding countries, the locations of Edinburgh, Bergen, and the Harding Oil and Gas field, the outline of Block 9 (the hydrocarbon exploration areas are divided into blocks, which are made up of the 1° longitudes and latitudes). It also includes the median line, which defines the economic sectors of the countries adjacent to the North Sea.
All data I use here are freely available online:
In [1]:
import shapefile
import numpy as np
from matplotlib import cm, rcParams
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
rcParams.update({'font.size': 16}) # Increase font-size
In [2]:
# Load the topo file to get header information
etopo1name = 'data/basemap/etopo1_bedrock.asc'
topo_file = open(etopo1name, 'r')
# Read header (number of columns and rows, cell-size, and lower left coordinates)
ncols = int(topo_file.readline().split()[1])
nrows = int(topo_file.readline().split()[1])
xllcorner = float(topo_file.readline().split()[1])
yllcorner = float(topo_file.readline().split()[1])
cellsize = float(topo_file.readline().split()[1])
topo_file.close()
# Read in topography as a whole, disregarding first five rows (header)
etopo = np.loadtxt(etopo1name, skiprows=5)
# Data resolution is quite high. I decrease the data resolution
# to decrease the size of the final figure
dres = 2
# Swap the rows
etopo[:nrows+1, :] = etopo[nrows+1::-1, :]
etopo = etopo[::dres, ::dres]
# Create longitude and latitude vectors for etopo
lons = np.arange(xllcorner, xllcorner+cellsize*ncols, cellsize)[::dres]
lats = np.arange(yllcorner, yllcorner+cellsize*nrows, cellsize)[::dres]
In [3]:
fig = plt.figure(figsize=(8, 6))
# Create basemap, 870 km east-west, 659 km north-south,
# intermediate resolution, Transverse Mercator projection,
# centred around lon/lat 1°/58.5°
m = Basemap(width=870000, height=659000,
resolution='i', projection='tmerc',
lon_0=1, lat_0=58.5)
# Draw coast line
m.drawcoastlines(color='k')
# Draw continents and lakes
m.fillcontinents(lake_color='b', color='none')
# Draw a think border around the whole map
m.drawmapboundary(linewidth=3)
# Convert etopo1 coordinates lon/lat in ° to x/y in m
# (From the basemap help: Calling a Basemap class instance with the arguments
# lon, lat will convert lon/lat (in degrees) to x/y map projection coordinates
# (in meters).)
rlons, rlats = m(*np.meshgrid(lons,lats))
# Draw etopo1, first for land and then for the ocean, with different colormaps
llevels = np.arange(-500,2251,100) # check etopo.ravel().max()
lcs = m.contourf(rlons, rlats, etopo, llevels, cmap=cm.terrain)
olevels = np.arange(-3500,1,100) # check etopo.ravel().min()
cso = m.contourf(rlons, rlats, etopo, olevels, cmap=cm.ocean)
# Draw parallels and meridians
m.drawparallels(np.arange(-56,63.,2.), color='.2', labels=[1,0,0,0])
m.drawparallels(np.arange(-55,63.,2.), color='.2', labels=[0,0,0,0])
m.drawmeridians(np.arange(-6.,12.,2.), color='.2', labels=[0,0,0,1])
m.drawmeridians(np.arange(-7.,12.,2.), color='.2', labels=[0,0,0,0])
# Draw Block 9 boundaries
m.plot([1, 2, 2, 1, 1], [59, 59, 60, 60, 59], 'b', linewidth=2, latlon=True)
plt.annotate('9', m(1.1, 59.7), color='b')
# Draw maritime boundaries
m.readshapefile('data/basemap/DECC_OFF_Median_Line', 'medline', linewidth=2)
# Add Harding, Edinburgh, Bergen
# 1. Convert coordinates
EDIx, EDIy = m(-3.188889, 55.953056)
BERx, BERy = m(5.33, 60.389444)
HARx, HARy = m(1.5, 59.29)
# 2. Plot symbol
plt.plot(HARx, HARy, mfc='r', mec='k', marker='s', markersize=10)
plt.plot(EDIx, EDIy, mfc='r', mec='k', marker='o', markersize=10)
plt.plot(BERx, BERy, mfc='r', mec='k', marker='o', markersize=10)
# 3. Plot name
plt.text(EDIx+50000, EDIy+10000,'Edinburgh', color='r')
plt.text(BERx-140000, BERy, 'Bergen', color='r')
plt.text(HARx-160000, HARy, 'Harding', color='r')
plt.show()