In the last two posts I show how to geo-reference a aerial images from various online services. I also commented that Google Earth Image was probably the best option from up-to-date image in small and/or dynamic areas. But what if you need bathymetry data? What if your image is cover by ice?
When all global datasets fails you can always resort to good old Nautical Charts. Below I show how to download, georeference, and plot a nautical chart using basemap.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
def make_map(llcrnrlon=None, urcrnrlon=None,
llcrnrlat=None, urcrnrlat=None,
img=None, figsize=(10, 10), resolution='c'):
m = Basemap(llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon,
llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat,
projection='merc', resolution=resolution,
lat_ts=(llcrnrlat + urcrnrlat) / 2.)
fig, ax = plt.subplots(figsize=figsize, facecolor='none')
m.ax = ax
if img:
image = plt.imread(img)
m.imshow(image, origin='upper', alpha=0.75)
meridians = np.linspace(llcrnrlon, urcrnrlon, 4)
parallels = np.linspace(llcrnrlat, urcrnrlat, 4)
kw = dict(linewidth=0)
m.drawparallels(parallels, labels=[1, 0, 0, 0], **kw)
m.drawmeridians(meridians, labels=[0, 0, 0, 1], **kw)
return fig, m
I will download a BSB chart from the
Brazilian navy website.
The chart files are zip, so we need to unzip, convert from BSB to PNG and
sub-sample them to make the figure lighter.
Below there is a function to do all that. (CAVEAT: Sometimes the Brazilian navy packs the files inside a sub-directory, sometimes not, sometimes there are more files than just the charts, sometimes not...)
import os
from zipfile import ZipFile
from urllib.request import urlretrieve
def get_chart(chart_number, outdir="."):
# Download.
fname = '%s/chart.zip' % outdir
uri = "https://www.mar.mil.br/dhn/chm/cartas/download/cartasbsb"
url = "%s/%s.zip" % (uri, number)
urlretrieve(url, fname)
# Extract.
with ZipFile(fname) as zfile:
BSB = zfile.namelist()[0]
KAP = zfile.namelist()[1]
zfile.extract(BSB, outdir)
zfile.extract(KAP, outdir)
PNG = '%s/%s.png' % (outdir, os.path.splitext(KAP)[0])
KAP = './data/%s' % KAP
if os.system('bsb2png %s %s' % (KAP, PNG)) == 0:
print("Conversion succeeded, file %s created." % PNG)
if os.system('convert %s -resize 10%% %s' % (PNG, PNG)) == 0:
print("File size reduced.")
return PNG
Here I download the chart number 25121 for the Admiralty bay in Antarctica.
number = 25121
outdir = './data'
PNG = get_chart(number, outdir=outdir)
The longitude and latitude corners are written in the image. If you download a different chart, just open the image with your favorite image viewer and change the values below.
llcrnrlon = -(58 + 40/60)
llcrnrlat = -(62 + 18/60)
urcrnrlon = -(58 + 12/60)
urcrnrlat = -(62 + 2/60)
Now the final cell! An geo-referenced nautical chart with an global inset showing the big picture.
fig, m = make_map(llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon,
llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat,
figsize=(8, 8), img=PNG, resolution='c')
m.fillcontinents()
m.drawcoastlines()
if True:
# Global inset map.
axin = inset_axes(m.ax, width=2, height=2, loc=4,
bbox_to_anchor=(0.85, 0.15),
bbox_transform=m.ax.figure.transFigure)
inmap = Basemap(projection='ortho',
lon_0=np.mean([llcrnrlon, urcrnrlon]),
lat_0=-90, ax=axin, anchor='NE')
inmap.drawcountries(color='white')
inmap.fillcontinents(color='gray')
bx, by = inmap(m.boundarylons, m.boundarylats)
xy = list(zip(bx, by))
mapboundary = Polygon(xy, edgecolor='k', linewidth=1, fill=False)
inmap.ax.add_patch(mapboundary)
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.