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

"""
Explanation: Plotting the location of the North Sea Harding Oil and Gas field
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:
<ul>
  <li> `etopo1_bedrock.asc`: You can download the etopo-data from the
    U.S. [National Geophysical Data Center (NGDC)](http://maps.ngdc.noaa.gov/viewers/wcs-client/).
    For this example, I downloaded a data set with the following parameters:</li>
  <ol>
    <li> ETOPO1 (bedrock)</li>
    <li> -8 West to 10 East, 64 North to 54 South</li>
    <li> ArcGIS ASCII Grid</li>
  </ol>
  <li> `DECC_OFF_Median_Line`: The median line is downloaded from the
    UK [Department of Energy & Climate Change (DECC)](https://www.gov.uk/oil-and-gas-offshore-maps-and-gis-shapefiles).
    The coordinates of Harding are from Well 9/23b-7, which are also available at the
    [DECC](https://www.gov.uk/oil-and-gas-wells).</li>
  <li> The coordinates of Edinburgh and Bergen are taken from Wikipedia.</li>
</ul>
End of explanation
"""


# 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]

"""
Explanation: Load data
End of explanation
"""


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()

"""
Explanation: Create the figure
End of explanation
"""


