<jupyter_start><jupyter_text>read in ECCO LLC data and regrid to uniform grid for saving as netcdf and plotting- first pip install stuff & restart kernal<jupyter_code>#!pip install git+https://github.com/rabernat/xmitgcm.git@fix-llcreader-klevels-bug
!pip install ecco_v4_py<jupyter_output><empty_output><jupyter_text>import libs<jupyter_code>import numpy as np
import cmocean
import xarray as xr
import ecco_v4_py as ecco
import xmitgcm.llcreader as llcreader
%matplotlib inline
import holoviews as hv
#from holoviews.operation.datashader import regrid
#hv.extension('bokeh')

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (15,10)

%matplotlib inline
import holoviews as hv
from holoviews.operation.datashader import regrid
hv.extension('bokeh')<jupyter_output><empty_output><jupyter_text>skip this unless you need to read in original data<jupyter_code>model = llcreader.ECCOPortalLLC2160Model()
ds = model.get_dataset(k_levels=[0], type='latlon')
ds
## read in data and get Theta
model = llcreader.ECCOPortalLLC2160Model()
ds_sst = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon')
ds_wnd = model.get_dataset(varnames=['oceTAUX','oceTAUY'], type='latlon')
ds_wnd<jupyter_output><empty_output><jupyter_text>regrid data onto new uniform grid using eccopy<jupyter_code>ecco_ds = ds_sst
new_grid_delta_lat = 360/8640
new_grid_delta_lon = 360/8640

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.Theta.isel(time=6200),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

da = xr.DataArray(field_nearest_1deg,name='Theta',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))
da.to_netcdf('./../data/ecco2011_11_19_latlon360Theta.nc')
da.coords['lon'] = np.mod(da['lon'], 360)
da = da.sortby(da.lon)
# output data
da.to_netcdf('./../data/ecco2011_11_19_latlonTheta.nc')<jupyter_output><empty_output><jupyter_text>switch from -180-180 to 0-360 to look at pacific ocean<jupyter_code>ecco_ds = ds_wnd
new_grid_delta_lat = 360/8640
new_grid_delta_lon = 360/8640

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.oceTAUX.isel(time=6200),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)
da = xr.DataArray(field_nearest_1deg,name='oceTAUX',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))
da.to_netcdf('./../data/ecco2011_11_19_latlon360taux.nc')
da.coords['lon'] = np.mod(da['lon'], 360)
da = da.sortby(da.lon)
# output data
da.to_netcdf('./../data/ecco2011_11_19_latlontaux.nc')

ecco_ds = ds_wnd
new_grid_delta_lat = 360/8640
new_grid_delta_lon = 360/8640

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.oceTAUY.isel(time=6200),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)
da = xr.DataArray(field_nearest_1deg,name='oceTAUY',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))
da.to_netcdf('./../data/ecco2011_11_19_latlon360tauy.nc')
da.coords['lon'] = np.mod(da['lon'], 360)
da = da.sortby(da.lon)
# output data
da.to_netcdf('./../data/ecco2011_11_19_latlontauy.nc')<jupyter_output><empty_output><jupyter_text>start from here if just plotting existing data<jupyter_code>import cartopy
import cartopy.crs as ccrs
import xarray as xr
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
da = xr.open_dataset('./../data/ecco2011_11_19_latlon.nc')
da1 = xr.open_dataset('./../data/ecco2011_11_19_latlontaux.nc')
da2 = xr.open_dataset('./../data/ecco2011_11_19_latlontauy.nc')
da['spd']= np.sqrt(da1.oceTAUX**2+da2.oceTAUY**2)
data = da['spd']
data = np.nan_to_num(data,0)
lowpass = gaussian_filter(data, 40)
gauss_highpass = data - lowpass
da['spd2'] = xr.DataArray(gauss_highpass,coords={'lat':da.lat,'lon':da.lon},dims=('lat','lon'))
da['spd3'] = xr.DataArray(lowpass,coords={'lat':da.lat,'lon':da.lon},dims=('lat','lon'))

#plt.rcParams['figure.figsize'] = (7,5)
plt.figure(figsize=(18,7))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = da.sel(lat=slice(-50,-35),lon=slice(10,90)).Theta.plot(ax=ax, 
                            cmap='turbo', 
                            vmin=2, 
                            vmax=20) # plot a colormap in transformed coordinates cmap='OrRd', 
plt.savefig('./../figures/theta_oneill_ortho_turbo.png')<jupyter_output><empty_output><jupyter_text>exp with wind below ehre<jupyter_code>lat1,lat2 = -46,-38
lon1,lon2 = 52,59
subset = da.sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
X, Y = np.meshgrid(subset.lon,subset.lat)
plt.figure(figsize=(15,7))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = subset.Theta.plot(ax=ax,cmap='turbo',vmin=3,vmax=18) # plot a colormap in transformed coordinates cmap='OrRd', 
q2 = subset.spd2.plot.contour(ax = ax, levels=np.arange(-.2,.2,.02))
#cp = ax.contourf(X, Y, subset.spd2)
lat1,lat2 = -46,-38
lon1,lon2 = 52,59
subset = da.sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
subset2 = subset.interp(lat=np.arange(subset.lat[0],subset.lat[-1],.25),lon=np.arange(subset.lon[0],subset.lon[-1],.25),method='nearest')
X, Y = np.meshgrid(subset2.lon,subset2.lat)
plt.figure(figsize=(15,7))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = subset2.Theta.plot(ax=ax,cmap='turbo',vmin=3,vmax=18) # plot a colormap in transformed coordinates cmap='OrRd', 
q2 = subset2.spd2.plot.contour(ax = ax, levels=np.arange(-.2,.2,.02))
#cp = ax.contourf(X, Y, subset.spd2)
lat1,lat2 = -46,-38
lon1,lon2 = 52,59
subset = da.sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
subset2 = subset.interp(lat=np.arange(subset.lat[0],subset.lat[-1],.5),lon=np.arange(subset.lon[0],subset.lon[-1],.5),method='nearest')
X, Y = np.meshgrid(subset2.lon,subset2.lat)
plt.figure(figsize=(15,7))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = subset2.Theta.plot(ax=ax,cmap='turbo',vmin=3,vmax=18) # plot a colormap in transformed coordinates cmap='OrRd', 
q2 = subset2.spd2.plot.contour(ax = ax, levels=np.arange(-.2,.2,.02))
#cp = ax.contourf(X, Y, subset.spd2)
lat1,lat2 = -46,-38
lon1,lon2 = 52,59
subset = da.sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
X, Y = np.meshgrid(subset.lon,subset.lat)
plt.figure(figsize=(15,7))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = subset.Theta.plot(ax=ax,cmap='turbo',vmin=3,vmax=18) # plot a colormap in transformed coordinates cmap='OrRd', 
q2 = subset.spd2.plot.contour(ax = ax, levels=np.arange(-.2,.2,.02))
#cp = ax.contourf(X, Y, subset.spd2)

lat1,lat2 = -60,-40
lon1,lon2 = 270,360
subset = da.sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
X, Y = np.meshgrid(subset.lon,subset.lat)
plt.figure(figsize=(9,7))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = subset.Theta.plot(ax=ax,cmap='turbo',vmin=0,vmax=20) # plot a colormap in transformed coordinates cmap='OrRd', 
q2 = subset.spd2.plot.contour(ax = ax, levels=np.arange(-.2,.2,.05))

# try an average

new_grid_delta_lat = 360/8640
new_grid_delta_lon = 360/8640

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

## read in data and get Theta
model = llcreader.ECCOPortalLLC2160Model()
ds_sst = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon')
ds_wnd = model.get_dataset(varnames=['oceTAUX','oceTAUY'], type='latlon')

ecco_ds = ds_sst
new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.Theta.isel(time=slice(6200,6230)),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

                                
da = xr.DataArray(field_nearest_1deg,name='Theta',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))

ecco_ds = ds_wnd
new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.oceTAUX.isel(time=slice(6200,6230)),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)
da2 = xr.DataArray(field_nearest_1deg,name='oceTAUX',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.oceTAUY.isel(time=slice(6200,6230)),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)
da3 = xr.DataArray(field_nearest_1deg,name='oceTAUY',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))




ecco_ds = ds_wnd
ecco_ds
plt.figure(figsize=(18,14))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = da.sel(lat=slice(-60,-28),lon=slice(10,90)).spd.plot(ax=ax, 
                            cmap='turbo', 
                            vmin=0, 
                            vmax=1) # plot a colormap in transformed coordinates cmap='OrRd', 
#try to lowpass filter data

plt.figure(figsize=(18,14))                    #set the figure size
ax = plt.subplot(1, 1, 1) #, projection=ortho)  #create the axis for plotting
q = da.sel(lat=slice(-60,-28),lon=slice(10,90)).spd2.plot(ax=ax, 
                            cmap='turbo', 
                            vmin=-.1, 
                            vmax=.1) # plot a colorma

ortho = ccrs.Orthographic(-170, 20)           # define target coordinate frame
geo = ccrs.PlateCarree()                     # define origin coordinate frame

plt.figure(figsize=(18,14))                    #set the figure size
ax = plt.subplot(1, 1, 1, projection=ortho)  #create the axis for plotting

q = da.Theta.plot(ax=ax, 
                            transform = geo, 
                            cmap='turbo', 
                            vmin=5, 
                            vmax=29) # plot a colormap in transformed coordinates cmap='OrRd', 

#global_extent = ax.get_extent(crs=ccrs.PlateCarree())
#gg = global_extent[:2] + (-60,-28)
#gg = (10,90)+gg[2:] 
#ax.set_extent(gg, crs=ccrs.PlateCarree())
#ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.COASTLINE)
plt.savefig('./../figures/theta_global2_ortho_turbo.png')
import cartopy
import cartopy.crs as ccrs
ortho = ccrs.Orthographic(-170, 20)           # define target coordinate frame
geo = ccrs.PlateCarree()                     # define origin coordinate frame

plt.figure(figsize=(18,14))                    #set the figure size
ax = plt.subplot(1, 1, 1, projection=ortho)  #create the axis for plotting

q = da.Theta.plot(ax=ax, 
                            transform = geo, 
                            cmap=cmocean.cm.thermal, 
                            vmin=5, 
                            vmax=29) # plot a colormap in transformed coordinates cmap='OrRd', 

ax.add_feature(cartopy.feature.COASTLINE)
plt.savefig('./../figures/theta_global2_ortho_thermal.png')<jupyter_output><empty_output><jupyter_text>you can also just plot data even if not on a regular grid<jupyter_code>plt.rcParams['figure.figsize'] = (20,10)
ds_sst.Theta.isel(k=0,time=6200,j=slice(2500,6480),i=slice(3000,8000)).plot(vmin=5,vmax=29,cmap='turbo')
plt.savefig('./../figures/theta_global.png')
plt.rcParams['figure.figsize'] = (20,10)
ds_sst.Theta.isel(k=0,time=6200,j=slice(2500,6480),i=slice(3000,8000)).plot(vmin=5,vmax=28,cmap='turbo')
plt.savefig('./../figures/theta_global2.png')
plt.rcParams['figure.figsize'] = (20,10)
ds_sst.Theta.isel(k=0,time=6200,j=slice(2500,6480),i=slice(3000,8000)).plot(vmin=5,vmax=28,cmap='coolwarm')
plt.savefig('./../figures/theta_global2_coolwarm.png')

plt.rcParams['figure.figsize'] = (20,10)
ds_sst.Theta.isel(k=0,time=6200,j=slice(2500,6480),i=slice(3000,8000)).plot(vmin=5,vmax=28,cmap=cmocean.cm.thermal)
plt.savefig('././figures/theta_global2_thermal.png')<jupyter_output><empty_output><jupyter_text>EXPLORE DATA<jupyter_code>ds_all = model.get_dataset(k_levels=[0], type='latlon')
ds_all
plt.rcParams['figure.figsize'] = (20,10)
ds_all.oceQnet.isel(time=6200,j=slice(2500,6480),i=slice(3000,8000)).plot(vmin=-200,vmax=400,cmap=cmocean.cm.thermal)
plt.savefig('./../figures/theta_global2_heatflux_thermal2.png')
#testing for when to look
dataset = hv.Dataset(ds_sst.Theta.isel(k=0,j=slice(3500,6100),i=slice(3000,8000),time=slice(1,18000,240)).astype('f4'))
hv_im = (dataset.to(hv.Image, ['i', 'j'], dynamic=True)
                .options(cmap='Magma', width=950, height=600, colorbar=True))

%output holomap='scrubber' fps=3
regrid(hv_im, precompute=True)<jupyter_output><empty_output><jupyter_text>make more data netcdf files read in data and get Theta<jupyter_code>model = llcreader.ECCOPortalLLC2160Model()
ds_sst = model.get_dataset(varnames=['oceQnet'], k_levels=[0], type='latlon')

ecco_ds = ds_sst
new_grid_delta_lat = 360/8640
new_grid_delta_lon = 360/8640

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.oceQnet.isel(time=6200),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

da = xr.DataArray(field_nearest_1deg,name='oceQnet',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))
# output data
da.to_netcdf('./data/ecco2011_11_19_latlon360oceQnet.nc')
da
# try to make 3D
## read in data and get Theta
model = llcreader.ECCOPortalLLC2160Model()
ds_sst = model.get_dataset(varnames=['Theta'], type='latlon')
ds_sst
ds_sst = ds_sst.isel(time=6200)
ds_sst
tem = ds_sst.isel(i=slice(4000,4200),j=slice(1000,1200))
tem
ecco_ds = ds_sst
new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.Theta,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

                                
da = xr.DataArray(field_nearest_1deg,name='Theta',coords={'lat':new_grid_lat[:,0],'lon':new_grid_lon[0,:]},dims=('lat','lon'))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1372439582519395, 'fraction_numerical': 0.05006418485237484, 'mean_word_length': 2.9782415630550623, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '252006', 'license': 'Apache-2.0', 'n_tokens_mistral': 7372, 'n_tokens_neox': 6716, 'n_words': 939, 'path': 'Butterfly/ecco_figure_sst_wnd.ipynb', 'repo_name': 'sevfour/cloud_science', 'size': '17917'}