Scott Lindell at MBL was working on aquaculture siting in Cape Cod Bay, and wanted to know when and where endangered Right Whales were found in shallow water.
In [26]:
import pandas as pd
import netCDF4
from mpl_toolkits.basemap import interp
Scott's data was in a simple spreadsheet output, so we saved to CSV and read it into Pandas
In [27]:
df = pd.read_csv('/usgs/data2/notebook/right_whale.csv',index_col='DATE',parse_dates=True)
In [28]:
df.head(4)
Out[28]:
To see what depths the whales were in, get a high-resolution bathymetry dataset for the Gulf of Maine. We will get just the bathymetric grid we need by determining the bounding box of the whale data, and then using that to extract the data via OPeNDAP. We access OPeNDAP data using NetCDF4-Python, which makes OPeNDAP datasets look like a local NetCDF file.
In [29]:
bathy_url='http://geoport.whoi.edu/thredds/dodsC/bathy/gom03_v1_0'
nc = netCDF4.Dataset(bathy_url).variables
lon=nc['lon'][:]
lat=nc['lat'][:]
bi=(lon>=df['Longitude'].min())&(lon<=df['Longitude'].max())
bj=(lat>=df['Latitude'].min())&(lat<=df['Latitude'].max())
z=nc['topo'][bj,bi]
lon=lon[bi]
lat=lat[bj]
print shape(z)
Get whale depths by interpolating from bathymetric grid
In [30]:
df['Depth']=interp(-z,lon,lat,df['Longitude'],df['Latitude'])
In [31]:
df.head()
Out[31]:
At what water depths are whales usually found in Cape Cod Bay?
In [32]:
hist(df['Depth']);
In [33]:
# create dataframe of shallow sightings
ishallow = where(df['Depth'] < 10.0)[0]
df_shallow = df.ix[ishallow]
df_shallow.head(10)
Out[33]:
In [34]:
figure(figsize=(12,12))
z = ma.masked_where(z>0,z)
subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
pcolormesh(lon,lat,z,vmax=0)
plot(df_shallow.Longitude,df_shallow.Latitude,'ko');
plot(df_shallow.Longitude,df_shallow.Latitude,'w+');
axis([-70.6, -70.0, 41.75, 42.06])
grid()
In [35]:
# determine stats by year (or month, quarter, etc)
#df_shallow_stats = df_shallow.resample('6M', how='mean') # monthly means
df_shallow_stats = df_shallow.resample('A-AUG', how='sum') # annual sum, ending in august
In [36]:
df_shallow_stats['SPECIES TOTAL']
Out[36]:
In [37]:
# any significant trend in time?
df_shallow_stats['SPECIES TOTAL'].plot(kind='bar',figsize=(12,5))
Out[37]:
So how many whales on average would we expect to see in shallow water in specific months? We can use the Pandas groupby method to group by month and take the mean
In [64]:
# average number of sightings for each month
df_monthly_sum = df_shallow.groupby(lambda x: x.month).mean()
In [67]:
df_monthly_sum.head(12)
Out[67]:
In [68]:
# plot 6 frame panel of sightings in each month
iframe=0
figure(figsize=(20,10))
grouped = df_shallow.groupby(lambda x: x.month)
for month,group in grouped:
iframe = iframe + 1
subplot(2,3,iframe,aspect=(1.0/cos(mean(lat)*pi/180.0)))
pcolormesh(lon,lat,z,vmax=0)
plot(group.Longitude,group.Latitude,'ko')
plot(group.Longitude,group.Latitude,'w+')
axis('tight')
#axis('off')
grid('on')
title('Month = %d' % month)
In [42]:
for month,group in grouped:
group.to_csv('/usgs/data2/notebook/%s.csv' % month)
In [43]:
# plot each month separately
iframe=0
for month,group in grouped:
iframe = iframe + 1
#subplot(2,3,iframe,aspect=(1.0/cos(mean(lat)*pi/180.0)))
figure(iframe,figsize=(12,12))
subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
pcolormesh(lon,lat,z,vmax=0)
plot(group.Longitude,group.Latitude,'ko')
plot(group.Longitude,group.Latitude,'w+')
axis([-70.6, -70.0, 41.75, 42.06])
#axis('tight')
#axis('off')
grid('on')
title('Month = %d' % month)
In [55]:
df_shallow_stats[['SPECIES TOTAL','Depth']].describe()
Out[55]:
In [45]:
In [46]:
df_shallow.index
Out[46]:
In [47]:
type(df_shallow.index)
Out[47]:
In [ ]: