In [2]:
%matplotlib inline
import warnings
import matplotlib.cbook
warnings.filterwarnings("ignore",category=matplotlib.cbook.mplDeprecation)
whales get entangled
Whales get sotted by planes, ships, and gliders.
https://whalesandmarinefauna.files.wordpress.com/2013/01/628x471_004.jpg
https://media.mnn.com/assets/images/2017/12/TangledNorthAtlanticRightWhale.jpg
spoter planes https://pbs.twimg.com/media/DdVa8ltX0AEo10V.jpg
https://assets.skiesmag.com/images/online/MaritimeProtector/1.jpg
https://assets.skiesmag.com/images/online/MaritimeProtector/4.jpg
http://www.releases.gov.nl.ca/releases/2013/fishaq/1216n04_pic.jpg
https://ssl.c.photoshelter.com/img-get2/I0000jNrkh7emj4Q/fit=1000x750/GSBCAN0401005.jpg
http://imgproc.airliners.net/photos/airliners/5/7/1/2711175.jpg?v=v40
https://twitter.com/FishOceansCAN/status/991054848538923008
Spoter vessels CCGS Corporal McLaren M.M.V. https://navaltoday.com/wp-content/uploads/2015/05/Canadian-Coast-Guard-Expands-Fleet.jpg
http://www.ccg-gcc.gc.ca/Vessel-Procurement/Gallery
Present the data (I think it is better to do this bit by bit... present each data set just before it is used)
WHale possitions
PLane surveys
Glider surveys
In [3]:
import pandas as pd
whales = pd.read_csv('whale_sightings.csv')
print(whales)
basic "map-style" plotting (non projected)
In [4]:
whales.plot(x='Longitude',y='Latitude',kind='scatter')
Out[4]:
Basemap
In [5]:
from mpl_toolkits.basemap import Basemap
# Create map
m = Basemap(projection='mill',
llcrnrlat=39.9,
urcrnrlat=48.3,
llcrnrlon=-69,
urcrnrlon=-54.7,
resolution='l')
x, y = m(whales['Longitude'].values,whales['Latitude'].values)
cs = m.scatter(x,y,s=20,marker='o',color='r')
m.fillcontinents()
Out[5]:
In [6]:
def make_basemap():
from mpl_toolkits.basemap import Basemap
# Create map
m = Basemap(projection='mill',
llcrnrlat=39.9,
urcrnrlat=48.3,
llcrnrlon=-69,
urcrnrlon=-54.7,
resolution='l')
m.fillcontinents()
return m
def plot_whales(whales, m):
x, y = m(whales['Longitude'].values,whales['Latitude'].values)
cs = m.scatter(x,y,s=40,marker='o',color='r',zorder=50)
return
In [7]:
make_basemap()
Out[7]:
In [8]:
make_basemap()
plot_whales(whales, m)
In [9]:
def make_basemap():
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
fig1 = plt.figure(figsize=(20,14))
# Create map
m = Basemap(projection='mill',
llcrnrlat=39.9,
urcrnrlat=48.3,
llcrnrlon=-69,
urcrnrlon=-54.7,
resolution='i')
m.drawcoastlines(color='#a6a6a6',linewidth=0.5,zorder=42)
m.fillcontinents(color='#e6e6e6',zorder=42)
m.drawmapboundary()
return m
In [10]:
m = make_basemap()
plot_whales(whales, m)
In [11]:
def add_maritimes_region(m):
import numpy as np
import shapefile
sf = shapefile.Reader("M:\\OCMD\\01_Data\\Management_Areas\\Administrative_Boundaries\\RegionalPolygon_Maritimes_OCMD_2013\\OriginalData\\MaritimesRegionPolygon_UpdatedSept2015_wgs84")
for shape in list(sf.iterShapes()):
npoints=len(shape.points) # total points
nparts = len(shape.parts) # total parts
if nparts == 1:
poly_lons = np.zeros((len(shape.points),1))
poly_lats = np.zeros((len(shape.points),1))
for ip in range(len(shape.points)):
poly_lons[ip] = shape.points[ip][0]
poly_lats[ip] = shape.points[ip][1]
plot_polygon(poly_lons, poly_lats)
else: # loop over parts of each shape, plot separately
for ip in range(nparts): # loop over parts, plot separately
i0=shape.parts[ip]
if ip < nparts-1:
i1 = shape.parts[ip+1]-1
else:
i1 = npoints
seg=shape.points[i0:i1+1]
poly_lons = np.zeros((len(seg),1))
poly_lats = np.zeros((len(seg),1))
for ip in range(len(seg)):
poly_lons[ip] = seg[ip][0]
poly_lats[ip] = seg[ip][1]
plot_polygon(poly_lons, poly_lats,m, edgecolor='#ff0000',linewidth=1.0,alpha=0.6,zorder=40)
return
def plot_polygon(poly_lons, poly_lats, m, edgecolor='#a6a6a6',linewidth=0.5,alpha=0.3,zorder=2):
import numpy as np
from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
poly_x, poly_y = m(poly_lons, poly_lats)
poly_xy = np.transpose(np.array((poly_x[:,0], poly_y[:,0])))
# Ad polygon
poly = Polygon(poly_xy,
closed=True,
edgecolor=edgecolor,
linewidth=linewidth,
alpha=alpha,
fill=False,
zorder=zorder)
plt.gca().add_patch(poly)
return
In [12]:
m = make_basemap()
plot_whales(whales, m)
add_maritimes_region(m)
In [ ]:
def add_NAFO_areas(m):
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import shapefile
sf = shapefile.Reader("M:\\OCMD\\01_Data\\Management_Areas\\Fisheries\\FishingZones_NEAtlantic_NAFO_2011\\ValueAdded\\NAFO_SubUnits_CanAtlantic")
for shape in list(sf.iterShapes()):
npoints=len(shape.points) # total points
nparts = len(shape.parts) # total parts
if nparts == 1:
poly_lons = np.zeros((len(shape.points),1))
poly_lats = np.zeros((len(shape.points),1))
for ip in range(len(shape.points)):
poly_lons[ip] = shape.points[ip][0]
poly_lats[ip] = shape.points[ip][1]
plot_polygon(poly_lons, poly_lats,m, edgecolor='#a6a6a6',linewidth=0.5,alpha=0.5,zorder=20)
else: # loop over parts of each shape, plot separately
for ip in range(nparts): # loop over parts, plot separately
i0=shape.parts[ip]
if ip < nparts-1:
i1 = shape.parts[ip+1]-1
else:
i1 = npoints
seg=shape.points[i0:i1+1]
poly_lons = np.zeros((len(seg),1))
poly_lats = np.zeros((len(seg),1))
for ip in range(len(seg)):
poly_lons[ip] = seg[ip][0]
poly_lats[ip] = seg[ip][1]
plot_polygon(poly_lons, poly_lats, m, edgecolor='#a6a6a6',linewidth=0.5,alpha=0.5,zorder=20)
# NAFO labels
nafo = pd.read_csv('M:\\OCMD\\02_Projects\\IOM\\NARW_ResourceMgmt\\NAFO_subunit_centroids.csv')
zones = pd.unique(nafo['UnitArea'].values)
for zone in zones:
zone_points = nafo[nafo['UnitArea'] == zone]
lat = zone_points['ddlat'].values[0]
lon = zone_points['ddlong'].values[0]
if lat > 39.9 and lat < 48.3 and lon > -69 and lon < -54.7:
plot_label(lon, lat, zone, m)
return
def plot_label(lon, lat, zone, m):
import matplotlib.pyplot as plt
x, y = m(lon, lat)
plt.text(x, y, zone, fontsize=9,color='#a6a6a6',zorder=35)
return
In [ ]:
m = make_basemap()
#plot_whales(whales,m)
add_maritimes_region(m)
add_NAFO_areas(m)
In [13]:
def add_bathymetry(m):
import matplotlib.pyplot as plt
import netCDF4
import urllib.request
import numpy as np
bathymetry_file = 'usgsCeSrtm30v6.nc'
minlat=39.9
maxlat=48.3
minlon=-69
maxlon=-54.7
isub = 1
base_url='http://coastwatch.pfeg.noaa.gov/erddap/griddap/usgsCeSrtm30v6.nc?'
query='topo[(%f):%d:(%f)][(%f):%d:(%f)]' % (maxlat,isub,minlat,minlon,isub,maxlon)
url = base_url+query
# store data in NetCDF file
urllib.request.urlretrieve(url, bathymetry_file)
# open NetCDF data in
nc = netCDF4.Dataset(bathymetry_file)
ncv = nc.variables
lon = ncv['longitude'][:]
lat = ncv['latitude'][:]
lons, lats = np.meshgrid(lon, lat)
topo = ncv['topo'][:, :]
TOPOmasked = np.ma.masked_where(topo>0,topo)
# For topo
x, y = m(lons, lats)
plt.pcolormesh(x,y,TOPOmasked,cmap=plt.get_cmap('Blues_r'),zorder=5,vmax=2000)
depth_levels_1 = np.linspace(topo.min(), -700, num=5)
depth_levels = np.append(depth_levels_1,np.linspace(-650, -50, num=15))
depth_levels = depth_levels.tolist()
cs = plt.contour(
x,
y,
topo,
depth_levels,
cmap=plt.get_cmap('Blues_r'),
linewidths=0.3,
linestyles='solid',
zorder=19)
return
Make a basemap with bathymetry, and features
In [14]:
m = make_basemap()
#plot_whales(whales,m)
#add_maritimes_region(m)
#add_NAFO_areas(m)
add_bathymetry(m)
Save basemap with pickle
OPen basemap, replot whales
In [15]:
#vms_csv = r'C:\Users\IbarraD\Data\VMS_DFO_Oracle\data_original\vms_autoDownloaded.csv'
vms_csv = 'vms_autoDownloaded.csv'
fishing = pd.read_csv(vms_csv)
print(fishing)
OPen basemap, replot VMS
In [16]:
import matplotlib.pyplot as plt
m = make_basemap()
plot_whales(whales,m)
#add_maritimes_region(m)
#add_NAFO_areas(m)
add_bathymetry(m)
x, y = m(fishing['LONGITUDE'].values,fishing['LATITUDE'].values)
plt.scatter(x,y,s=1,marker='o',color='y', zorder=10)
Out[16]:
In [42]:
def make_heatmap(info, m, data):
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
print('gridder ---------------------------------------------')
x = np.arange(info['minlon'], info['maxlon'], info['bin_size'], dtype=np.float64)
y = np.arange(info['minlat'], info['maxlat'], info['bin_size'], dtype=np.float64)
info['bin_number'] = (np.ceil((info['maxlon'] - info['minlon'])/info['bin_size']),
np.ceil((info['maxlat'] - info['minlat'])/info['bin_size']))
x,y = m(data['LONGITUDE'].values, data['LATITUDE'].values)
# Project pings to grid
#H, xedges, yedges = np.histogram2d(x,y,bins=info['bin_number'],
# range=[[0, info['bin_number'][0]],
# [0, info['bin_number'][1]]])
H, xedges, yedges = np.histogram2d(x,y,bins=500)
# Rotate and flip H...
H = np.rot90(H)
H = np.flipud(H)
# Mask zeros
Hmasked = np.ma.masked_where(H==0,H)
# Log H for better display
Hmasked = np.log10(Hmasked)
# Make map
#plt.pcolormesh(xedges[0:-1],yedges[0:-1],H,zorder=1000)
# plt.show()
cs = m.pcolor(xedges,yedges,Hmasked,cmap=plt.get_cmap('inferno_r'), zorder=1000)
# plt.show()
return
In [43]:
#from gridder import make_heatmap
info = {}
info['minlat'] = 39.9
info['maxlat'] = 48.3
info['minlon'] = -69
info['maxlon'] = -54.7
info['bin_size'] = 0.01 # Degrees
m = make_basemap()
plot_whales(whales,m)
#add_maritimes_region(m)
#add_NAFO_areas(m)
add_bathymetry(m)
make_heatmap(info, m, fishing)
In [ ]:
In [ ]:
In [ ]:
In [19]:
import importlib
importlib.reload(make_heatmap)
Filter by species
Filter by speed
Install ship_mapper
Use ship_mapper to make fishing density heatmap
Supperimpose whale locations
In [ ]: