Location of high salinity points

Created by Ivan Lima on Wed Jul 22 2020 16:06:51 -0400

In [1]:
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import os, datetime, warnings, cmocean
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from cesm_utils import fixtime, da2ma
from tqdm.notebook import trange, tqdm
warnings.filterwarnings('ignore')
print('Last updated on {}'.format(datetime.datetime.now().ctime()))
%matplotlib inline
plt.rcParams['figure.dpi'] = 110

Last updated on Wed Jul 22 16:10:29 2020
In [2]:
infile = '/media/ivan/WD/Data/Postproc/1DPOP/cesm.1dpop.07.abs.max_2D.nc'
ds_max = xr.open_dataset(infile, decode_times=False,
                         drop_variables=['EVAP_F','ROFF_F','IFRAC','HMXL','TEMP','UVEL','VVEL','ULAT','ULONG','z_w','dz','dzw'])
ds_max = fixtime(ds_max, 2051)

stacked = ds_max.SALT.stack(z=('nlat','nlon'))
lonmax, latmax = [], []
imax, jmax = [], []
ntime = ds_max.dims['time']
for n in trange(1,11):
    maxind = np.argsort(stacked.fillna(-999), axis=1)[:,-n]
    maxpos = stacked['z'][maxind]
    lonmax = lonmax + [maxpos.TLONG[t].item() for t in range(ntime)]
    latmax = latmax + [maxpos.TLAT[t].item() for t in range(ntime)]
    imax = imax + [maxpos.values[t][0] for t in range(ntime)]
    jmax = jmax + [maxpos.values[t][1] for t in range(ntime)]



















In [3]:



    


fig, ax = plt.subplots(subplot_kw={'projection':ccrs.Robinson(central_longitude=210)})
ax.set_global()
ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
ax.coastlines(linewidth=0.5)
for lon, lat in sorted(set([p for p in zip(lonmax,latmax)])):
    print('lat = {:.2f}N, lon = {:.2f}E'.format(lat, lon))
    _ = ax.plot(lon, lat, 'C3x', transform=ccrs.PlateCarree(), alpha=0.75)

_ = ax.set_title('position of salinity maxima')
# fig.savefig('salt_max_pos.png', dpi=200)



    















    






lat = 37.04N, lon = 11.58E
lat = 36.57N, lon = 11.62E
lat = 36.10N, lon = 11.65E
lat = 36.68N, lon = 13.84E
lat = 36.20N, lon = 13.88E
lat = 40.68N, lon = 24.65E
lat = 40.75N, lon = 25.76E
lat = 31.44N, lon = 29.86E
lat = 31.52N, lon = 32.10E
lat = 36.27N, lon = 35.16E
lat = 29.51N, lon = 49.06E
lat = 28.98N, lon = 49.09E
lat = 28.45N, lon = 49.11E
lat = 27.92N, lon = 49.13E
lat = 29.55N, lon = 50.19E
lat = 29.02N, lon = 50.21E
lat = 28.48N, lon = 50.23E
lat = 27.95N, lon = 50.25E
lat = 29.05N, lon = 51.34E
lat = 28.52N, lon = 51.36E
lat = 31.33N, lon = 122.44E
lat = 60.13N, lon = 314.85E
lat = 60.09N, lon = 316.72E
lat = 60.47N, lon = 316.75E










    


























In [4]:



    


for i, j in sorted(set([p for p in zip(imax,jmax)])):
    print('({},{:3d}) lat = {:.2f}N, lon = {:.2f}E'.format(i, j, ds_max.TLAT[i,j].item(), ds_max.TLONG[i,j].item()))



    















    






(271, 79) lat = 27.92N, lon = 49.13E
(271, 80) lat = 27.95N, lon = 50.25E
(272, 79) lat = 28.45N, lon = 49.11E
(272, 80) lat = 28.48N, lon = 50.23E
(272, 81) lat = 28.52N, lon = 51.36E
(273, 79) lat = 28.98N, lon = 49.09E
(273, 80) lat = 29.02N, lon = 50.21E
(273, 81) lat = 29.05N, lon = 51.34E
(274, 79) lat = 29.51N, lon = 49.06E
(274, 80) lat = 29.55N, lon = 50.19E
(274,144) lat = 31.33N, lon = 122.44E
(279, 62) lat = 31.44N, lon = 29.86E
(279, 64) lat = 31.52N, lon = 32.10E
(288, 67) lat = 36.27N, lon = 35.16E
(290, 46) lat = 36.10N, lon = 11.65E
(290, 48) lat = 36.20N, lon = 13.88E
(291, 46) lat = 36.57N, lon = 11.62E
(291, 48) lat = 36.68N, lon = 13.84E
(292, 46) lat = 37.04N, lon = 11.58E
(298, 58) lat = 40.68N, lon = 24.65E
(298, 59) lat = 40.75N, lon = 25.76E
(353,314) lat = 60.13N, lon = 314.85E
(353,316) lat = 60.09N, lon = 316.72E
(354,316) lat = 60.47N, lon = 316.75E



















In [5]:



    


minlon, maxlon = 300, 330
minlat, maxlat =  52, 65
lonc = (minlon + maxlon)/2.
latc = (minlat + maxlat)/2.
proj = ccrs.EquidistantConic(central_longitude=lonc,central_latitude=latc)
# proj = ccrs.Mercator(central_longitude=lonc)

tlat, tlon = ds_max.TLAT, ds_max.TLONG
tlat, tlon = np.c_[tlat,tlat[:,0]], np.c_[tlon,tlon[:,0]]

fig, ax = plt.subplots(subplot_kw={'projection':proj}, figsize=(9,6.4))
ax.set_extent([minlon,maxlon,minlat,maxlat])
ax.coastlines(linewidth=0.5)
# ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
gl = ax.gridlines(xlocs=np.arange(-180,180,10),ylocs=np.arange(0,90,5),draw_labels=True)
gl.xlabels_top, gl.ylabels_right = False, False
pm = ax.pcolormesh(tlon, tlat, ds_max.SALT[-1,:,:], transform=ccrs.PlateCarree(), cmap=cmocean.cm.haline,
                   vmin=25, vmax=54.25)
cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.075)
cb.set_label(ds_max.SALT.units)
_ = ax.set_title('salinity on %d-%02d-%02d'%(ds_max.time[-1]['time.year'],
                                                     ds_max.time[-1]['time.month'],
                                                     ds_max.time[-1]['time.day']))
for lon, lat in sorted(set([p for p in zip(lonmax,latmax)])):
    _ = ax.plot(lon, lat, 'C3o', transform=ccrs.PlateCarree(), alpha=0.75)



    















    


























In [6]:



    


minlon, maxlon = -7, 37
minlat, maxlat =  28, 45
lonc = (minlon + maxlon)/2.
latc = (minlat + maxlat)/2.
proj = ccrs.EquidistantConic(central_longitude=lonc,central_latitude=latc)

fig, ax = plt.subplots(subplot_kw={'projection':proj}, figsize=(9,6.4))
ax.set_extent([minlon,maxlon,minlat,maxlat])
ax.coastlines(linewidth=0.5)
# ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
gl = ax.gridlines(xlocs=np.arange(-180,180,10),ylocs=np.arange(0,90,5),draw_labels=True)
gl.xlabels_top, gl.ylabels_right = False, False
pm = ax.pcolormesh(tlon, tlat, ds_max.SALT[-1,:,:], transform=ccrs.PlateCarree(), cmap=cmocean.cm.haline,
                   vmin=36, vmax=43)
cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.075)
cb.set_label(ds_max.SALT.units)
_ = ax.set_title('salinity on %d-%02d-%02d'%(ds_max.time[-1]['time.year'],
                                                     ds_max.time[-1]['time.month'],
                                                     ds_max.time[-1]['time.day']))
for lon, lat in sorted(set([p for p in zip(lonmax,latmax)])):
    _ = ax.plot(lon, lat, 'C3o', transform=ccrs.PlateCarree(), alpha=0.75)



    















    


























In [10]:



    


minlon, maxlon = 40, 63
minlat, maxlat =  9, 32
lonc = (minlon + maxlon)/2.
latc = (minlat + maxlat)/2.
proj = ccrs.EquidistantConic(central_longitude=lonc,central_latitude=latc)

fig, ax = plt.subplots(subplot_kw={'projection':proj}, figsize=(9,6.4))
ax.set_extent([minlon,maxlon,minlat,maxlat])
ax.coastlines(linewidth=0.5)
# ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
gl = ax.gridlines(xlocs=np.arange(-180,180,10),ylocs=np.arange(0,90,5),draw_labels=True)
gl.xlabels_top, gl.ylabels_right = False, False
pm = ax.pcolormesh(tlon, tlat, ds_max.SALT[-1,:,:], transform=ccrs.PlateCarree(), cmap=cmocean.cm.haline,
                   vmin=35, vmax=43)
cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.075)
cb.set_label(ds_max.SALT.units)
_ = ax.set_title('salinity on %d-%02d-%02d'%(ds_max.time[-1]['time.year'],
                                                     ds_max.time[-1]['time.month'],
                                                     ds_max.time[-1]['time.day']))
for lon, lat in sorted(set([p for p in zip(lonmax,latmax)])):
    _ = ax.plot(lon, lat, 'C3o', transform=ccrs.PlateCarree(), alpha=0.75)