In [1]:
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
import bq
import time
In [2]:
def get_area(lat):
lat_degree = 69 #miles
# Convert latitude and longitude to
# spherical coordinates in radians.
degrees_to_radians = math.pi/180.0
# phi = 90 - latitude
phi = (lat+.05)*degrees_to_radians #plus 0.5 to get the middle
lon_degree = math.cos(phi)*69 #miles
# return 69*69*2.6
return lat_degree*lon_degree* 2.58999 #miles to square km
In [3]:
client = bq.Client.Get()
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 [4]:
# query to get the density of vessels
q = '''
SELECT integer(lat*10) lat_bin, integer(lon*10) lon_bin,
sum(1/positions) vessels
FROM [scratch_david_gapanalysis.orbcomm_2015_noduplicates_st_stats]
where
max_lat - min_lat <5
AND (max_lon - min_lon < 10
OR first_lon > 170
OR first_lon < -170)
and lat > -90 and lat < 90 and lat != 0 and lon > -180 and lon < 180
AND mmsi IN (select mmsi from
[scratch_david_gapanalysis.good_mmsi_2015])
//and tagblock_type = 'terrestrial'
group by lat_bin, lon_bin
'''
vessel_query = Query(q)
In [12]:
# set all pixels to zero that have a value, and the rest to 1
vessel_days = np.ones(shape=(1800,3600))
# zero anywhere there is data
for row in vessel_query:
lat = int(row[0])
lon = int(row[1])
if lat<900 and lat>-900 and lon>-1800 and lon<1800:
lat_index = lat+900
lon_index = lon+1800
vessel_days[lat_index][lon_index] = 0
In [19]:
firstlat = 90
lastlat = -90
firstlon = -180
lastlon = 180
scale = .1
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()
m = Basemap(llcrnrlat=lastlat, urcrnrlat=firstlat,
llcrnrlon=firstlon, urcrnrlon=lastlon, lat_ts=0, projection='mill',resolution="h")
m.drawmapboundary()
# m.drawcoastlines(linewidth=.2)
#m.fillcontinents('#555555')#, lake_color, ax, zorder, alpha)
x = np.linspace(-180, 180, 360*10 )
y = np.linspace(lastlat, firstlat, (firstlat-lastlat)*10)
x, y = np.meshgrid(x, y)
converted_x, converted_y = m(x, y)
from matplotlib import colors,colorbar
# norm = colors.LogNorm(vmin=1, vmax=1000)
norm = colors.Normalize(vmin=0, vmax=1)
m3 = 1
m.pcolormesh(converted_x, converted_y, vessel_days, norm=norm, vmin=0, vmax=1)
#plt.pcolormesh(x, y, z, cmap='RdBu', vmin=z_min, vmax=z_max)
t = "Where there be no boats, 2015"
plt.title(t)
ax = fig.add_axes([0.15, 0.1, 0.4, 0.02]) #x coordinate ,
norm = colors.Normalize(vmin=0, vmax=1)
# norm = colors.LogNorm(vmin=1, vmax=1000)
lvls = np.logspace(0,3,7)
cb = colorbar.ColorbarBase(ax,norm = norm, orientation='horizontal', ticks=lvls)
cb.ax.set_xticklabels(["<1" ,int(m3**.5), m3, int(m3**1.5), m3*m3,int(m3**2.5), str(int(m3**3))+"+"], fontsize=10)
cb.set_label('vessels - yes or no',labelpad=-40, y=0.45)
plt.savefig("vessel_density_2015_zeros.png",bbox_inches='tight',dpi=300,transparent=True,pad_inches=0)
plt.rcParams["figure.figsize"] = [12,8]
plt.show()
#plt.clf()
In [20]:
v = np.zeros((1800))
for i in range(1800):
v[i] = vessel_days[i][3030:3060].sum()
plt.scatter([i for i in range(1800)], v)
In [15]:
plt.clf()
plt.imshow(np.flipud(vessel_days), interpolation='nearest')#, norm=norm)
plt.savefig("vessel_density_2015_zeros_imshow.png",bbox_inches='tight',dpi=300,transparent=True,pad_inches=0)
In [ ]: