Fishing Effort in the White Shark Cafe


In [3]:
import bq  
import time
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib
from matplotlib import colors,colorbar
import matplotlib
%matplotlib inline
import csv 
import math
from math import radians, cos, sin, asin, sqrt
from scipy import stats
import matplotlib.dates as mdates

client = bq.Client.Get()

In [4]:
# create a bounding box:
max_lat = 27
min_lat = 20
max_lon = -131
min_lon = -135

In [5]:
def Query(q):
    t0 = time.time()
    answer = client.ReadTableRows(client.Query(q)['configuration']['query']['destinationTable'])
    print 'Query time: ' + str(time.time() - t0) + ' seconds.'
    return answer

In [6]:
q = '''
select count(distinct mmsi) 
FROM
  [tilesets.pipeline_2015_08_24_08_19_01]
WHERE
  latitude > '''+str(min_lat)+'''
  AND latitude <'''+str(max_lat)+'''
  AND longitude > '''+str(min_lon)+'''
  AND longitude < '''+str(max_lon)+'''
  AND weight >=.5'''

number_of_mmsi = Query(q)


Waiting on bqjob_r36f9187a7778e44f_000001545e85dc2e_1 ... (48s) Current status: DONE   
Query time: 50.740391016 seconds.

In [7]:
print "Number of unique MMSI:",int(number_of_mmsi[0][0])


Number of unique MMSI: 6

6 MMSI were fishing in this region from January 2014 to July 2015


In [12]:
q = '''
select a.mmsi mmsi,
a.number number,
if(b.country is null, "invalid mmsi",b.country) country,
b.continent continent
from
(select mmsi, 
integer(if(length(string(mmsi))= 9,LEFT(STRING(mmsi),3), '0')) code,
count(*) number from
(SELECT
  mmsi,
  date(timestamp) date,
  first(latitude) lat,
  first(longitude) lon
FROM
  [tilesets.pipeline_2015_08_24_08_19_01]
WHERE
  latitude > '''+str(min_lat)+'''
  AND latitude <'''+str(max_lat)+'''
  AND longitude > '''+str(min_lon)+'''
  AND longitude < '''+str(max_lon)+'''
  AND weight >=.5
group by mmsi, date) 
group by mmsi, code 
) 
a
left join [scratch_roan.country_code] b
on a.code = b.code 
order by number desc'''

days_by_mmsi = Query(q)


Waiting on bqjob_r3d088056c9986cff_000001545e88bc7b_3 ... (22s) Current status: DONE   
Query time: 24.9715480804 seconds.

The List of MMSI that were fishing are below


In [13]:
for r in days_by_mmsi:
    if "Taiwan" in r[2]:
        r[2] = "Taiwan"
    print "mmsi:", r[0],"  fishing days:", r[1], "  country: ", r[2].split(" (")[0]


mmsi: 367338850   fishing days: 3   country:  United States of America
mmsi: 367512170   fishing days: 2   country:  United States of America
mmsi: 416001473   fishing days: 1   country:  Taiwan
mmsi: 9103923   fishing days: 1   country:  invalid mmsi
mmsi: 416004507   fishing days: 1   country:  Taiwan
mmsi: 367721000   fishing days: 1   country:  United States of America

In [14]:
# Now Group by Country
q = '''
SELECT 
country,
sum(number) number
from
(select a.mmsi mmsi,
a.number number,
if(b.country is null, "invalid mmsi",b.country) country,
b.continent continent
from
(select mmsi, 
integer(if(length(string(mmsi))= 9,LEFT(STRING(mmsi),3), '0')) code,
count(*) number from
(SELECT
  mmsi,
  date(timestamp) date,
  first(latitude) lat,
  first(longitude) lon
FROM
  [tilesets.pipeline_2015_08_24_08_19_01]
WHERE
  latitude > '''+str(min_lat)+'''
  AND latitude <'''+str(max_lat)+'''
  AND longitude > '''+str(min_lon)+'''
  AND longitude < '''+str(max_lon)+'''
  AND weight >=.5
group by mmsi, date) 
group by mmsi, code 
) 
a
left join [scratch_roan.country_code] b
on a.code = b.code)
group by country
order by number desc'''

country_groups = Query(q)


Waiting on bqjob_r60f3eb3c52c83615_000001545e8a1f6f_4 ... (22s) Current status: DONE   
Query time: 25.0199449062 seconds.

In [15]:
print "days\tcountry"
for c in country_groups:
    print c[1],'\t', c[0].split(" (")[0]


days	country
6 	United States of America
2 	Taiwan
1 	invalid mmsi

In [16]:
objects = [c[0].split(" (")[0] for c in country_groups]
y_pos = np.arange(len(objects))
performance = [c[1] for c in country_groups]
 
plt.figure(figsize=(12,4))
plt.bar(y_pos-4, performance, align='center', alpha=0.5,)
plt.xticks(y_pos-4+.2, objects, fontsize=12) # rotation=45
locs, labels = plt.xticks()
plt.setp(labels, rotation=45,h`a='right')
plt.ylabel('Usage',fontsize=12)
plt.title('Fishing Days by Country in Mascarene Plateau, Jan 2014 to July 2015',fontsize=15)

ax = plt.axes()
ax.xaxis.set_ticks_position('none') 
plt.savefig("fishing_effort_cafe_by_Country.png",bbox_inches='tight',dpi=150,transparent=True,pad_inches=.1)
plt.show()


Map the Fishing Effort in This Region


In [17]:
# Now Group by Country
q = '''
SELECT
lat,
lon,
count(*) fishing_days
from
(SELECT
  mmsi,
  date(timestamp) date,
  integer(first(latitude)*4) lat,
  integer(first(longitude)*4) lon
FROM
  [tilesets.pipeline_2015_08_24_08_19_01]
WHERE
  latitude > '''+str(min_lat)+'''
  AND latitude <'''+str(max_lat)+'''
  AND longitude > '''+str(min_lon)+'''
  AND longitude < '''+str(max_lon)+'''
  AND weight >=.5
group by mmsi, date) 
group by lat, lon
'''

fishing_grid = Query(q)


Waiting on bqjob_r2e337a2f67a2d5b8_000001545e8a83b7_5 ... (34s) Current status: DONE   
Query time: 36.0829241276 seconds.

In [18]:
cellsize = .25
one_over_cellsize = 4

num_lats = (max_lat-min_lat)*one_over_cellsize+1
num_lons = (max_lon-min_lon)*one_over_cellsize+1


grid = np.zeros(shape=(num_lats,num_lons))

for row in fishing_grid:
    lat = int(row[0])
    lon = int(row[1])
    lat_index = lat-min_lat*one_over_cellsize
    lon_index = lon-min_lon*one_over_cellsize
    grid[lat_index][lon_index] = int(row[2])

In [19]:
plt.rcParams["figure.figsize"] = [8,9]

firstlat = max_lat
lastlat =  min_lat
firstlon = min_lon
lastlon = max_lon
scale = cellsize

numlats = int((firstlat-lastlat)/scale+.5)
numlons = int((lastlon-firstlon)/scale+.5)
    
lat_boxes = np.linspace(lastlat,firstlat,num=numlats,endpoint=False)
lon_boxes = np.linspace(firstlon,lastlon,num=numlons,endpoint=False)

fig = plt.figure()
extra = 10
m = Basemap(llcrnrlat=lastlat-extra, urcrnrlat=firstlat+extra,
          llcrnrlon=firstlon-extra, urcrnrlon=lastlon+extra, lat_ts=0, projection='mill',resolution="h")

m.drawmapboundary()#fill_color='#111111')
m.drawcoastlines(linewidth=.2)
m.fillcontinents('#cccccc',lake_color='#cccccc')#, lake_color, ax, zorder, alpha)

x = np.linspace(firstlon, lastlon, -(firstlon-lastlon)*one_over_cellsize+1)
y = np.linspace(lastlat, firstlat, (firstlat-lastlat)*one_over_cellsize+1)
x, y = np.meshgrid(x, y)
converted_x, converted_y = m(x, y)
from matplotlib import colors,colorbar

maximum = grid.max()
minimum = 1

norm = colors.LogNorm(vmin=minimum, vmax=maximum)
# norm = colors.Normalize(vmin=0, vmax=1000)

m.pcolormesh(converted_x, converted_y, grid, norm=norm, vmin=minimum, vmax=maximum, cmap = plt.get_cmap('viridis'))

t = "Fishing Days Jan 2014 to July 2015"
plt.title(t, color = "#000000", fontsize=18)

ax = fig.add_axes([0.2, 0.1, 0.65, 0.02]) #x coordinate , 
norm = colors.LogNorm(vmin=minimum, vmax=maximum)
# norm = colors.Normalize(vmin=0, vmax=1000)
lvls = np.logspace(np.log10(minimum),np.log10(maximum),num=8)
cb = colorbar.ColorbarBase(ax,norm = norm, orientation='horizontal', ticks=lvls, cmap = plt.get_cmap('viridis'))

#cb.ax.set_xticklabels(["0" ,round(m3**.5,1), m3, round(m3**1.5,1), m3*m3,round(m3**2.5,1), str(round(m3**3,1))+"+"], fontsize=10)
cb.ax.set_xticklabels([int(i) for i in lvls], fontsize=10, color = "#000000")
cb.set_label('Fishing Days',labelpad=-40, y=0.45, color = "#000000")
plt.savefig("fishing_effort_cafe.png",bbox_inches='tight',dpi=300,transparent=True,pad_inches=.1)
plt.show()



In [ ]: