In [37]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
In [1]:
import matplotlib.pyplot as plt
%matplotlib inline
import cartopy
import cartopy.crs as ccrs
import numpy as np
import xarray as xr
import math
In [2]:
amm7ds = xr.open_dataset('data/amm7_subdomains_v2.nc')
cp2_lv1_fname = 'data/CP2_lv1_v2.shp' # One region
cp2_lv2_fname = 'data/CP2_lv2_v2.shp' # UK shelf anf off shelf
cp2_lv3_fname = 'data/CP2_new_edit_v2.shp'# The complete 9 regions
In [3]:
def plot_shapefile(shapefile,extent_lonlat=[-16,6,46,65]):
plt.figure(figsize=(4,4))
crs_lambea = ccrs.LambertAzimuthalEqualArea(central_latitude=55.0)
crs_lonlat = ccrs.PlateCarree()
ax = plt.axes(projection=crs_lambea)
xmin, ymin = crs_lambea.transform_point(extent_lonlat[0], extent_lonlat[2], crs_lonlat)
xmax, ymax = crs_lambea.transform_point(extent_lonlat[1], extent_lonlat[3], crs_lonlat)
extent_lambea = (xmin, xmax, ymin, ymax)
ax.set_extent(extent_lambea,crs=crs_lambea)
shape_feature = ShapelyFeature(shapefile.geometries(),
ccrs.PlateCarree(), edgecolor='black')
ax.add_feature(shape_feature)#, facecolor='blue')
ax.coastlines(resolution='10m')
ax.gridlines()
In [47]:
# Plot facet-like all DataArrays in Dataset
def plot_ds(ds):
nvars = len(ds.data_vars)
nrows = math.ceil(math.sqrt(nvars))
ncols = math.ceil(nvars/nrows)
fig,axs = plt.subplots(nrows=nrows,ncols=ncols,figsize=(15,15),squeeze=False)
for ind, varn in enumerate(ds.data_vars):
data = ds[varn].values.squeeze()
col = ind % ncols
row = int((ind-col)/ncols)
#print('plot_ds',nrows,ncols,row,col,ind)
ax = axs[row,col].imshow(np.flipud(data))
axs[row,col].set_title(varn)
In [5]:
#**********************************************
def make_mask(polygon,arrayX,arrayY):
#**********************************************
[Nrows,Nlines] = arrayX.shape
mask = np.ma.zeros([Nrows,Nlines],dtype=bool)
start=time.time()
for row in range(Nrows):
for col in range(Nlines):
x = arrayX[row,col]
y = arrayY[row,col]
p=shply.Point(x,y)
if p.within(polygon):
mask[row,col] = True
stop=time.time()
elapsed= stop-start
elapseds="%.2f"%elapsed
print('make_mask: ' + elapseds + ' s')
return mask
In [40]:
# takes shpfile name and return list of raster masks
def shapefile2masks(shp_fname):
shpread = shp.Reader(shp_fname)
shapes = shpread.shapes()
#shapefile has no useful records describing the polygon
try:
shprecords = shpread.records()
#print('records:',shprecords)
shprecs=[]
for rec in shprecords:
# record 0 is OBJECTID
shprecs=rec.append(rec[0])
#print('shprecs',shprecs)
except:
print('ERROR: no records in shapefile')
masks = []
for shape in shapes:
#print(shape)
points = shape.points
poly = shply.Polygon(points)
mask = make_mask(poly,long,latg)
masks.append(mask)
return masks
In [6]:
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature
In [7]:
shapefile = shpreader.Reader(cp2_lv1_fname)
plot_shapefile(shapefile,extent_lonlat=[-16,6,46,65])
In [8]:
shapefile = shpreader.Reader(cp2_lv2_fname)
plot_shapefile(shapefile,extent_lonlat=[-16,6,46,65])
In [9]:
shapefile = shpreader.Reader(cp2_lv3_fname)
plot_shapefile(shapefile,extent_lonlat=[-16,6,46,65])
In [10]:
import shapefile as shp
import shapely.geometry as shply
import time
In [11]:
long,latg = np.meshgrid(amm7ds['lon'], amm7ds['lat'])
In [41]:
masks1 = shapefile2masks(cp2_lv1_fname)
masks2 = shapefile2masks(cp2_lv2_fname)
masks3 = shapefile2masks(cp2_lv3_fname)
In [48]:
# Get coordinates from AMM7 domain file
lat = amm7ds['lat']
lon = amm7ds['lon']
llcoords = {'lat':lat,'lon':lon}
lldims = ('lat','lon')
In [50]:
# Create Dataset with all masks and lat lon dimensions and save to netCDF
da0 = xr.DataArray(data=masks1[0],coords=llcoords,dims=lldims)
cp2_lv1_ds = xr.Dataset({'CP2_lv1':da0})
# Save to NetCDF file
#cp2_lv2_ds.to_netcdf('data/amm7_CP2_lv2.nc')
plot_ds(cp2_lv1_ds)
In [51]:
# Create Dataset with all masks and lat lon dimensions and save to netCDF
print(len(masks2))
#plt.imshow(np.flipud(masks1[0]))
da0 = xr.DataArray(data=masks2[0],coords=llcoords,dims=lldims)
da1 = xr.DataArray(data=masks2[1],coords=llcoords,dims=lldims)
da2 = xr.DataArray(data=masks2[2],coords=llcoords,dims=lldims)
cp2_lv2_ds = xr.Dataset({'CP2_lv2_0':da0,'CP2_lv2_1':da1,'CP2_lv2_2':da2})
# Save to NetCDF file
#cp2_lv2_ds.to_netcdf('data/amm7_CP2_lv2.nc')
plot_ds(cp2_lv2_ds)
In [52]:
fileds = xr.open_dataset('data/amm7_CP2_lv2.nc')
id0 = fileds['CP2_lv2_0']
id1 = fileds['CP2_lv2_1']
id2 = fileds['CP2_lv2_2']
In [53]:
# Combine polygons to test for overlap or gaps
fig = plt.figure( figsize=(10,10))
(id1+id0+id2).plot()
Out[53]:
In [58]:
# create Dataset with all masks and lat lon dimensions
# and save to netCDF
print(len(masks3))
da0 = xr.DataArray(data=masks3[0],coords=llcoords,dims=lldims)
da1 = xr.DataArray(data=masks3[1],coords=llcoords,dims=lldims)
da2 = xr.DataArray(data=masks3[2],coords=llcoords,dims=lldims)
da3 = xr.DataArray(data=masks3[3],coords=llcoords,dims=lldims)
da4 = xr.DataArray(data=masks3[4],coords=llcoords,dims=lldims)
da5 = xr.DataArray(data=masks3[5],coords=llcoords,dims=lldims)
da6 = xr.DataArray(data=masks3[6],coords=llcoords,dims=lldims)
da7 = xr.DataArray(data=masks3[7],coords=llcoords,dims=lldims)
da8 = xr.DataArray(data=masks3[8],coords=llcoords,dims=lldims)
da9 = xr.DataArray(data=masks3[9],coords=llcoords,dims=lldims)
cp2_lv3_ds = xr.Dataset({'Western Channel & Celtic Sea':da0,
'Eastern Channel':da1,
'Faroe-Shetland Channel':da2,
'Irish Sea':da3,'Minches & Western Scotland':da4,
'Northern North Sea':da5,
'Atlantic NW Approaches':da6,
'Rockall Trough & Bank':da7,
'Southern North Sea':da8,
'Scottish Continental Shelf':da9})
plot_ds(cp2_lv3_ds)
# Save to file
#cp2_lv3_ds.to_netcdf('data/amm7_CP2_lv3.nc')
In [59]:
fig = plt.figure( figsize=(10,10))
(da1+da0+da2+da3+da4+da5+da6+da7+da8+da9).plot()
Out[59]:
In [16]:
fileds = xr.open_dataset('data/amm7_CP2_lv3.nc')
fileds
#id0 = fileds['CP2_lv2_0']
#id1 = fileds['CP2_lv2_1']
#id2 = fileds['CP2_lv2_2']
Out[16]: