First of all, load the necessary libraries. These are the ones we discussed previously:
numpy
matplotlib
cartopy
xarray
In [1]:
importosimportnumpyasnpimportxarrayasxrimportcartopyimportcartopy.crsasccrsimportcartopy.featureascfeaturefromcartopy.mpl.gridlinerimportLONGITUDE_FORMATTER,LATITUDE_FORMATTERcartopy.config['data_dir']=os.getenv('CARTOPY_DIR',cartopy.config.get('data_dir'))importcmoceanfrommatplotlibimportpyplotasplt#%config InlineBackend.figure_format = 'retina'plt.ion()# To trigger the interactive inline modeimportwarningswarnings.filterwarnings("ignore",category=RuntimeWarning)
We will use the open_mfdataset function from xArray to open multiple netCDF files into a single xarray Dataset.
We will query load the GBR4km dataset from the AIMS server, so let's first define the base URL:
In [2]:
# For the bio datasetbase_url="http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/GBR4_H2p0_B3p1_Cq3b_Dhnd/daily-monthly/EREEFS_AIMS-CSIRO_GBR4_H2p0_B3p1_Cq3b_Dhnd_bgc_daily-monthly-"# For the hydro datasetbase_url2="http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-"
For the sake of the demonstration, we will only use 1 month:
In [4]:
month_st=1# Starting month month_ed=1# Ending month year=2018# Year# Based on the server the file naming convention biofiles=[f"{base_url}{year}-{month:02}.nc"formonthinrange(month_st,month_ed+1)]hydrofiles=[f"{base_url2}{year}-{month:02}.nc"formonthinrange(month_st,month_ed+1)]
Regridding of daily input data (from eReefs AIMS-CSIRO GBR4 BioGeoChemical 3.1 subset) from curvilinear (per input data) to rectilinear via inverse weighted distance from up to 4 closest cells.
ems_version :
v1.1.1 rev(6244M)
history :
Tue Oct 8 15:38:27 2019: ncatted -a positive,botz,o,char,up -a missing_value,botz,o,double,99. -a outside,botz,o,double,-9999. gbr4_bgc_all_simple_2018-01.nc
2020-08-20T23:45:30+10:00: vendor: AIMS; processing: None summaries
2020-08-21T23:07:30+10:00: vendor: AIMS; processing: None summaries
eReefs 4 km grid. SOURCE Catchments with 2019 condition from Dec 1, 2010 to June,30, 2018, Empirical SOURCE with 2019 condition, Jul 1, 2018 to April 30, 2019. More details of naming protocol at: eReefs.info.
:::{seealso}
Check out this link to see how to use Xarray’s convenient matplotlib-backed plotting interface to visualize your datasets!
:::
In [25]:
# Figure sizesize=(9,10)# Color from cmoceancolor=cmocean.cm.balance# Defining the figurefig=plt.figure(figsize=size,facecolor='w',edgecolor='k')# Axes with Cartopy projectionax=plt.axes(projection=ccrs.PlateCarree())# and extentax.set_extent([142.4,157,-7,-28.6],ccrs.PlateCarree())# Plotting using Matplotlib # We plot the PH at the surface at the final recorded time intervalcf=ds_bio.PH.isel(time=-1,k=-1).plot(transform=ccrs.PlateCarree(),cmap=color,vmin=7.975,vmax=8.125,add_colorbar=False)# Color barcbar=fig.colorbar(cf,ax=ax,fraction=0.027,pad=0.045,orientation="horizontal")cbar.set_label(ds_bio.PH.long_name+' '+ds_bio.PH.units,rotation=0,labelpad=5,fontsize=10)cbar.ax.tick_params(labelsize=8)# Titleplt.title('zc = '+str(ds_bio.zc.values.item(-1))+' m at '+str(ds_bio.coords['time'].values[-1])[:10],fontsize=11)# Plot lat/lon grid gl=ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,linewidth=0.1,color='k',alpha=1,linestyle='--')gl.top_labels=Falsegl.right_labels=Falsegl.xformatter=LONGITUDE_FORMATTERgl.yformatter=LATITUDE_FORMATTERgl.xlabel_style={'size':8}gl.ylabel_style={'size':8}# Add map features with Cartopy ax.add_feature(cfeature.NaturalEarthFeature('physical','land','10m',edgecolor='face',facecolor='lightgray'))ax.coastlines(linewidth=1)# Site Davies Reefax.scatter(reef_lon,reef_lat,c='deeppink',s=50,edgecolors='k',linewidth=1,transform=ccrs.PlateCarree())plt.tight_layout()plt.show()fig.clear()plt.close(fig)plt.clf()
To reduce the Dataset size we will clip the spatial extent based on longitudinal and latitudinal values.
This is easely done using the sel function with the slice method.
In [15]:
min_lon=146# lower left longitudemin_lat=-21# lower left latitudemax_lon=150# upper right longitudemax_lat=-16# upper right latitude# Defining the boundarieslon_bnds=[min_lon,max_lon]lat_bnds=[min_lat,max_lat]# Performing the reductionds_bio_clip=ds_bio.sel(latitude=slice(*lat_bnds),longitude=slice(*lon_bnds))ds_hydro_clip=ds_hydro.sel(latitude=slice(*lat_bnds),longitude=slice(*lon_bnds))
Let's plot the clipped region using the same approach as above:
In [18]:
# Figure sizesize=(9,10)# Color from cmoceancolor=cmocean.cm.balance# Defining the figurefig=plt.figure(figsize=size,facecolor='w',edgecolor='k')# Axes with Cartopy projectionax=plt.axes(projection=ccrs.PlateCarree())# and extentax.set_extent([min_lon,max_lon,min_lat,max_lat],ccrs.PlateCarree())# Plotting using Matplotlib # We plot the PH at the surface at the final recorded time intervalcf=ds_bio_clip.PH.isel(time=-1,k=-1).plot(transform=ccrs.PlateCarree(),cmap=color,vmin=7.98,vmax=8.08,add_colorbar=False)# Color barcbar=fig.colorbar(cf,ax=ax,fraction=0.027,pad=0.045,orientation="horizontal")cbar.set_label(ds_bio_clip.PH.long_name+' '+ds_bio_clip.PH.units,rotation=0,labelpad=5,fontsize=10)cbar.ax.tick_params(labelsize=8)# Titleplt.title('zc = '+str(ds_bio_clip.zc.values.item(-1))+' m at '+str(ds_bio_clip.coords['time'].values[-1])[:10],fontsize=11)# Plot lat/lon grid gl=ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,linewidth=0.1,color='k',alpha=1,linestyle='--')gl.top_labels=Falsegl.right_labels=Falsegl.xformatter=LONGITUDE_FORMATTERgl.yformatter=LATITUDE_FORMATTERgl.xlabel_style={'size':8}gl.ylabel_style={'size':8}# Add map features with Cartopy ax.add_feature(cfeature.NaturalEarthFeature('physical','land','10m',edgecolor='face',facecolor='lightgray'))ax.coastlines(linewidth=1)# Site Davies Reefax.scatter(reef_lon,reef_lat,c='deeppink',s=70,edgecolors='k',linewidth=1,transform=ccrs.PlateCarree())plt.tight_layout()plt.show()fig.clear()plt.close(fig)plt.clf()
To make the quiver plot, we will first resample the dataset and take one point every 7 times. It allows to have less arrows on the map making it more readable.
In [20]:
# Figure sizesize=(7,8)# Time step to plottimevar=0# z-coordinate position (here the top one)zcvar=-1# Color from cmoceancolor=cmocean.cm.speed# Defining the figurefig=plt.figure(figsize=size,facecolor='w',edgecolor='k')# Axes with Cartopy projectionax=plt.axes(projection=ccrs.PlateCarree())# and extentax.set_extent([min_lon,max_lon,min_lat,max_lat],ccrs.PlateCarree())# Plotting using Matplotlib the mean currentcf=ds_hydro_clip.mean_cur.isel(time=timevar,k=zcvar).plot(transform=ccrs.PlateCarree(),cmap=color,vmin=0.1,vmax=0.7,add_colorbar=False)# Resampling using the slice methodresample=ds_hydro_clip.isel(time=timevar,k=zcvar,longitude=slice(None,None,7),latitude=slice(None,None,7))# Defining the quiver plotquiver=resample.plot.quiver(x='longitude',y='latitude',u='u',v='v',transform=ccrs.PlateCarree(),scale=8)# Vector options declarationveclenght=0.5maxstr='%3.1f m/s'%veclenghtplt.quiverkey(quiver,0.1,0.1,veclenght,maxstr,labelpos='S',coordinates='axes').set_zorder(11)# Color barcbar=fig.colorbar(cf,ax=ax,fraction=0.027,pad=0.045,orientation="horizontal")cbar.set_label(ds_hydro_clip.mean_cur.long_name+' '+ds_hydro_clip.mean_cur.units,rotation=0,labelpad=5,fontsize=10)cbar.ax.tick_params(labelsize=8)# Titleplt.title('zc = '+str(ds_hydro_clip.mean_cur.zc.values.item(zcvar))+' m at '+str(ds_hydro_clip.mean_cur.coords['time'].values[timevar])[:10],fontsize=11)# Plot lat/lon grid gl=ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,linewidth=0.1,color='k',alpha=1,linestyle='--')gl.top_labels=Falsegl.right_labels=Falsegl.xformatter=LONGITUDE_FORMATTERgl.yformatter=LATITUDE_FORMATTERgl.xlabel_style={'size':8}gl.ylabel_style={'size':8}# Add map features with Cartopy ax.add_feature(cfeature.NaturalEarthFeature('physical','land','10m',edgecolor='face',facecolor='lightgray'))ax.coastlines(linewidth=1)# Site Davies Reefax.scatter(reef_lon,reef_lat,c='deeppink',s=50,edgecolors='k',linewidth=1,transform=ccrs.PlateCarree()).set_zorder(11)plt.tight_layout()plt.show()fig.clear()plt.close(fig)plt.clf()
To define the cross-section, we chose to select (sel) several latitudes or longitudes and compute the mean along the opposite dimension:
In [22]:
# Section along the longitudelonsec=ds_bio_clip.sel(latitude=slice(*[-19.,-18.5])).mean(dim=('latitude'),skipna=True)# Section along the latitudelatsec=ds_bio_clip.sel(longitude=slice(*[147.2,148.])).mean(dim=('longitude'),skipna=True)
Let's visualise the selected mean cross-sections:
In [23]:
# Figure sizesize=(12,4)# Time step to plottimevar=-1# Color from cmoceancolor=cmocean.cm.speed# Defining the figurefig=plt.figure(figsize=size,facecolor='w',edgecolor='k')ax=plt.axes()# Plotting using Matplotlib the temperature over the z-coordinatesec=lonsec.temp.isel(time=timevar).plot(y="zc",add_colorbar=False)# Color barcbar=fig.colorbar(sec,ax=ax,fraction=0.027,pad=0.045)cbar.set_label(ds_bio.temp.long_name+' '+ds_bio.temp.units,rotation=90,labelpad=5,fontsize=10)cbar.ax.tick_params(labelsize=8)# Titleplt.title(str(lonsec.temp.coords['time'].values[timevar])[:10],fontsize=11)# Site Davies Reef visualised as a lineax.plot([reef_lon,reef_lon],[-1000,16],c='deeppink',linewidth=3)plt.tight_layout()plt.show()fig.clear()plt.close(fig)plt.clf()
<Figure size 432x288 with 0 Axes>
In [24]:
# Figure sizesize=(12,4)# Time step to plottimevar=-1# Color from cmoceancolor=cmocean.cm.speed# Defining the figurefig=plt.figure(figsize=size,facecolor='w',edgecolor='k')ax=plt.axes()# Plotting using Matplotlib the temperature over the z-coordinatesec=latsec.temp.isel(time=timevar).plot(y="zc",add_colorbar=False)# Color barcbar=fig.colorbar(sec,ax=ax,fraction=0.027,pad=0.045)cbar.set_label(ds_bio.temp.long_name+' '+ds_bio.temp.units,rotation=90,labelpad=5,fontsize=10)cbar.ax.tick_params(labelsize=8)# Titleplt.title(str(latsec.temp.coords['time'].values[timevar])[:10],fontsize=11)# Site Davies Reef visualised as a lineax.plot([reef_lat,reef_lat],[-1000,16],c='deeppink',linewidth=3)plt.tight_layout()plt.show()fig.clear()plt.close(fig)plt.clf()