In [2]:
%matplotlib inline
import warnings
import matplotlib.cbook
warnings.filterwarnings("ignore",category=matplotlib.cbook.mplDeprecation)

Whale sightings



In [3]:
import pandas as pd

whales = pd.read_csv('whale_sightings.csv')

print(whales)


         Date  Latitude  Longitude  Number of right whales
0  2018-05-31  47.59990 -64.098750                       1
1  2018-05-27  48.15185 -63.381437                       1
2  2018-05-24  44.58650 -63.058000                       1
3  2018-05-24  44.63517 -62.678800                       1

basic "map-style" plotting (non projected)



In [4]:
whales.plot(x='Longitude',y='Latitude',kind='scatter')


Out[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x1653a47e710>

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]:
[<matplotlib.patches.Polygon at 0x1653abec6d8>,
 <matplotlib.patches.Polygon at 0x1653abec9b0>,
 <matplotlib.patches.Polygon at 0x1653abecc50>,
 <matplotlib.patches.Polygon at 0x1653abecf28>]

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]:
<mpl_toolkits.basemap.Basemap at 0x1653a7999e8>

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)


---------------------------------------------------------------------------
ShapefileException                        Traceback (most recent call last)
<ipython-input-12-474f4e077033> in <module>()
      1 m = make_basemap()
      2 plot_whales(whales, m)
----> 3 add_maritimes_region(m)

<ipython-input-11-fa754e90616e> in add_maritimes_region(m)
      2     import numpy as np
      3     import shapefile
----> 4     sf = shapefile.Reader("M:\\OCMD\\01_Data\\Management_Areas\\Administrative_Boundaries\\RegionalPolygon_Maritimes_OCMD_2013\\OriginalData\\MaritimesRegionPolygon_UpdatedSept2015_wgs84")
      5     for shape in list(sf.iterShapes()):
      6         npoints=len(shape.points) # total points

~\Miniconda3workshop\lib\site-packages\shapefile.py in __init__(self, *args, **kwargs)
    235         if len(args) > 0:
    236             if is_string(args[0]):
--> 237                 self.load(args[0])
    238                 return
    239         if "shp" in kwargs.keys():

~\Miniconda3workshop\lib\site-packages\shapefile.py in load(self, shapefile)
    289                 pass
    290             if not (self.shp and self.dbf):
--> 291                 raise ShapefileException("Unable to open %s.dbf or %s.shp." % (shapeName, shapeName) )
    292         if self.shp:
    293             self.__shpHeader()

ShapefileException: Unable to open M:\OCMD\01_Data\Management_Areas\Administrative_Boundaries\RegionalPolygon_Maritimes_OCMD_2013\OriginalData\MaritimesRegionPolygon_UpdatedSept2015_wgs84.dbf or M:\OCMD\01_Data\Management_Areas\Administrative_Boundaries\RegionalPolygon_Maritimes_OCMD_2013\OriginalData\MaritimesRegionPolygon_UpdatedSept2015_wgs84.shp.

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



VMS



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)


       VR_NUMBER   LATITUDE  LONGITUDE    POSITION_UTC_DATE
0         105912  42.083000 -66.652000  2018-05-25 06:00:01
1         106604  42.130820 -66.917840  2018-05-25 06:00:01
2         107314  43.441650 -65.208410  2018-05-25 06:00:01
3         101019  43.744230 -65.317730  2018-05-25 06:00:01
4         100247  44.477891 -63.600112  2018-05-25 06:00:01
5         107918  44.478235 -63.599082  2018-05-25 06:00:01
6         103677  44.632472 -65.747166  2018-05-25 06:00:01
7         104994  44.729118 -63.008567  2018-05-25 06:00:01
8           7488  44.736580 -62.830900  2018-05-25 06:00:01
9         106311  44.907302 -65.867158  2018-05-25 06:00:01
10        103394  44.943523 -66.920558  2018-05-25 06:00:01
11        106991  44.949290 -65.265530  2018-05-25 06:00:01
12        107553  45.011210 -58.887130  2018-05-25 06:00:01
13        108065  45.053987 -66.798592  2018-05-25 06:00:01
14        103923  45.122180 -66.351250  2018-05-25 06:00:01
15        107204  45.231490 -60.181800  2018-05-25 06:00:01
16        106050  45.274228 -66.068259  2018-05-25 06:00:01
17        102680  45.338970 -60.997130  2018-05-25 06:00:01
18        108358  45.339031 -60.995494  2018-05-25 06:00:01
19        138113  45.339288 -60.997811  2018-05-25 06:00:01
20        108292  45.580558 -64.952632  2018-05-25 06:00:01
21        108029  45.582410 -60.740260  2018-05-25 06:00:01
22        104917  45.717170 -60.245190  2018-05-25 06:00:01
23        161024  45.899505 -59.976769  2018-05-25 06:00:01
24        106683  45.920440 -59.967890  2018-05-25 06:00:01
25        160602  46.006364 -59.840642  2018-05-25 06:00:01
26        158986  46.873420 -62.974880  2018-05-25 06:00:02
27         18437  45.338600 -60.995500  2018-05-25 06:00:15
28        104735  44.476740 -63.602280  2018-05-25 06:00:32
29        107181  44.482680 -63.594540  2018-05-25 06:00:32
...          ...        ...        ...                  ...
56354     153378  46.996164 -62.534180  2018-06-01 05:30:00
56355      12258  48.258670 -63.847730  2018-06-01 05:30:05
56356     139941  47.889180 -60.797340  2018-06-01 05:30:40
56357     103094  47.302890 -63.320260  2018-06-01 05:31:01
56358     153752  46.467000 -55.124000  2018-06-01 05:35:00
56359     153378  46.996850 -62.534695  2018-06-01 05:35:00
56360      12258  48.261970 -63.857830  2018-06-01 05:35:05
56361     153378  46.997451 -62.535296  2018-06-01 05:40:00
56362      12258  48.265280 -63.867930  2018-06-01 05:40:05
56363     151824  46.917500 -55.391000  2018-06-01 05:43:00
56364     133060  46.899003 -55.371094  2018-06-01 05:45:00
56365     153378  46.998052 -62.535897  2018-06-01 05:45:00
56366      12258  48.268470 -63.878030  2018-06-01 05:45:05
56367     139941  47.909160 -60.830460  2018-06-01 05:45:36
56368     103094  47.280660 -63.284970  2018-06-01 05:46:01
56369     153378  46.998653 -62.536412  2018-06-01 05:50:00
56370      12258  48.272280 -63.889030  2018-06-01 05:50:05
56371     153378  46.999254 -62.537013  2018-06-01 05:55:00
56372      12258  48.275180 -63.898130  2018-06-01 05:55:05
56373     100989  44.374894 -64.308558  2018-06-01 06:00:00
56374     131668  44.481300 -63.582370  2018-06-01 06:00:00
56375       5706  45.000000 -65.000000  2018-06-01 06:00:00
56376     155586  46.307910 -55.369960  2018-06-01 06:00:00
56377     102861  46.374290 -55.251230  2018-06-01 06:00:00
56378      84067  46.412700 -55.651810  2018-06-01 06:00:00
56379     137936  46.446570 -56.843030  2018-06-01 06:00:00
56380     101174  46.916700 -55.390130  2018-06-01 06:00:00
56381     153378  46.999854 -62.537614  2018-06-01 06:00:00
56382       5666  47.074441 -55.829001  2018-06-01 06:00:00
56383      94141  47.463400 -58.939490  2018-06-01 06:00:00

[56384 rows x 4 columns]

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]:
<matplotlib.collections.PathCollection at 0x1653e57e518>

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)


gridder ---------------------------------------------
C:\Users\cerc-user\Miniconda3workshop\lib\site-packages\ipykernel_launcher.py:32: RuntimeWarning: divide by zero encountered in log10

In [ ]:


In [ ]:


In [ ]:


In [19]:
import importlib
importlib.reload(make_heatmap)


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-19-1ea036724dbd> in <module>()
      1 import importlib
----> 2 importlib.reload(make_heatmap)

~\Miniconda3workshop\lib\importlib\__init__.py in reload(module)
    137     """
    138     if not module or not isinstance(module, types.ModuleType):
--> 139         raise TypeError("reload() argument must be a module")
    140     try:
    141         name = module.__spec__.name

TypeError: reload() argument must be a module

Filter by species



Filter by speed



Install ship_mapper



Use ship_mapper to make fishing density heatmap



Supperimpose whale locations



In [ ]: