# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # read in ECCO LLC data and regrid to uniform grid for saving as netcdf and plotting
#
# - first pip install stuff & restart kernal

# #!pip install git+https://github.com/rabernat/xmitgcm.git@fix-llcreader-klevels-bug
# !pip install ecco_v4_py

# import libs

# +
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')
# -

model = llcreader.ECCOPortalLLC2160Model()
ds = model.get_dataset(k_levels=[0], type='latlon')
ds

ds

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

# ## regrid data onto new uniform grid using eccopy

ds_sst.Theta.isel(time=-1,i=slice(5000,7000),j=slice(5000,6000)).plot()

# +
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=-1),\
                                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/ecco2013_04_22_latlon360Theta.nc')
da.coords['lon'] = np.mod(da['lon'], 360)
da = da.sortby(da.lon)
# output data
da.to_netcdf('./../data/ecco2013_04_22_latlonTheta.nc')
# -

# ## switch from -180-180 to 0-360 to look at pacific ocean

# ## start from here if just plotting existing data

import cartopy
import cartopy.crs as ccrs
import xarray as xr
import matplotlib.pyplot as plt
#da = xr.open_dataset('./../data/ecco2011_11_19_latlon.nc')
da = xr.open_dataset('./../data/ecco2013_04_22_latlonTheta.nc')

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)).Theta.plot(ax=ax, 
                            cmap='turbo', 
                            vmin=-1, 
                            vmax=22) # plot a colormap in transformed coordinates cmap='OrRd', 
plt.savefig('./../figures/theta_oneill_ortho_turbo.png')

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(45,63),lon=slice(180,242)).Theta.plot(ax=ax, 
                            cmap='turbo', 
                            vmin=2, 
                            vmax=9) # plot a colormap in transformed coordinates cmap='OrRd', 
plt.savefig('./../figures/theta_goa_ortho_turbo.png')

plt.figure(figsize=(18,7))                    #set the figure size
ax = plt.subplot(1, 1, 1, projection=ccrs.Orthographic(180, 10))  #create the axis for plotting
q = da.Theta.plot(ax=ax, 
                            cmap='turbo', 
                            transform=ccrs.PlateCarree(),
                            vmin=2, 
                            vmax=9) # plot a colormap in transformed coordinates cmap='OrRd', 
ax.add_feature(cartopy.feature.COASTLINE)
global_extent = ax.get_extent(crs=ccrs.PlateCarree())
gg = global_extent[:2] + (45,62)
gg = (180,242)+gg[2:] 
ax.set_extent(gg, crs=ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND)
plt.savefig('./../figures/theta_goa_ortho_turbo2.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='turbo', 
                            vmin=2, 
                            vmax=9) # plot a colormap in transformed coordinates cmap='OrRd', 

ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.COASTLINE)
plt.savefig('./../figures/theta_goa_ortho_turbo2.png')
# -

# ## you can also just plot data even if not on a regular grid

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')
# -

# ### EXPLORE DATA

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)
# -

# # make more data netcdf files
# ## read in data and get Theta
#

# +
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



StarcoderdataJupyterScriptsDedupFiltered