zachsa
commented
2 years ago Here is the Python script (from @mattcarr03) to create the example output:
- https://drive.google.com/file/d/1tqECbitcGe7Fh-Ibr9YyKPPkk-AMxm9x/view?usp=sharing
- https://drive.google.com/file/d/1b5dPq8GiTg-ENiCquo4cU5Wk6eDErpM_/view?usp=sharing
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 25 18:07:52 2022
@author: matthew
"""
#! /usr/bin/env python
#
# a script to animate croco surface currents with parcels particles
#
from matplotlib import pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
from netCDF4 import Dataset
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from datetime import timedelta, date, datetime
import sys
#sys.path.insert(1, '/home/osboxes/2001_DEA/models/parcels/croco_test')
#import croco_tools_visualisation as croco
import xarray as xr
from shapely.geometry import Point, Polygon
import netCDF4 as nc
# %%
def v2rho_3d(var_v):
[T,M,Lp]=var_v.shape
var_rho=np.zeros((T,M+1,Lp))
var_rho[:,1:M-1,:]=0.5*np.squeeze([var_v[:,0:M-2,:]+var_v[:,1:M-1,:]])
var_rho[:,0,:]=var_rho[:,1,:]
var_rho[:,M,:]=var_rho[:,M-1,:]
return var_rho
def u2rho_3d(var_u):
[T,Mp,L]=var_u.shape
var_rho=np.zeros((T,Mp,L+1))
var_rho[:,:,1:L-1]=0.5*np.squeeze([var_u[:,:,0:L-2]+var_u[:,:,1:L-1]])
var_rho[:,:,0]=var_rho[:,:,1]
var_rho[:,:,L]=var_rho[:,:,L-1]
return var_rho
def remove_child(lon_parent,lat_parent,data_parent,lon_child,lat_child):
# assign missing values for parent data inside child domain
# only works for regular grids
# get the lon lat coordinates along the boundary of the model (actually 1 grid cell in to ensure some overlap between parent and child grids)
lon_bry=lon_child[:,1].tolist()+lon_child[-2,:].tolist()+lon_child[-2:1:-1,-2].tolist()+lon_child[1,-2:1:-1].tolist()
lat_bry=lat_child[:,1].tolist()+lat_child[-2,:].tolist()+lat_child[-2:1:-1,-2].tolist()+lat_child[1,-2:1:-1].tolist()
#lon_bry=lon_child[:,0].tolist()+lon_child[-1,:].tolist()+lon_child[-1:0:-1,-1].tolist()+lon_child[0,-1:0:-1].tolist()
#lat_bry=lat_child[:,0].tolist()+lat_child[-1,:].tolist()+lat_child[-1:0:-1,-2].tolist()+lat_child[0,-1:0:-1].tolist()
bry = [(i,j) for i,j in zip(lon_bry,lat_bry)]
poly = Polygon(bry)
ny,nx=np.shape(lat_parent)
for j in range(ny):
for i in range(nx):
p = Point(lon_parent[j,i], lat_parent[j,i])
if p.within(poly):
data_parent[j,i]=np.nan
return data_parent
# %% define file names
path = '/home/matthew/confs/False_bay_forecast/croco/forecast/'
#child_path = '/home/matthew/confs/False_bay_forecast/croco/scratch/'
child_path = '/home/matthew/Downloads/'
# define file names
grdfile='grd.nc'
#avgchildfile='avg_20220225.nc.1'
avgchildfile = 'false-bay-forecast-20221101.nc'
#Loading reference dataset to find surface level index
ds_child = xr.open_dataset(child_path+avgchildfile)
# %% Loading variables from dataset
#Child
lon_child = ds_child.lon_rho.values
lat_child = ds_child.lat_rho.values
U_child = ds_child.u[:,19,:,:].values
V_child = ds_child.v[:,19,:,:].values
U_child = u2rho_3d(U_child)
V_child = v2rho_3d(V_child)
# %% Loading temperature
SST_child = ds_child.temperature[:,39,:,:].values
SST_child[np.where(SST_child==0)]=np.nan
# %% Loading time
# get the time
nc_id_avg = Dataset(child_path+avgchildfile, 'r')
time = nc_id_avg.variables['time'][:]
# time origin (needs to be hard coded here as it is not written to the croco output file)
date_ref=datetime(2000,1,1,0,0,0)
#generate a list of dates
dates=[]
for t in time:
date_now=date_ref + timedelta(seconds=np.float64(t))
dates.append(date_now)
# %% UV animation figure
skip = 2
# # setup figure object
fig = plt.figure(figsize=(12,8))
#ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
ax = plt.axes(projection=ccrs.PlateCarree())
#ax.set_extent([25.55,26.55,-34.1,-33.65])
ax.add_feature(cfeature.LAND, zorder=0, edgecolor='black')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle=':')
gl.left_labels = False
gl.top_labels = False
# add the vectors
U = U_child[0,:,:]
V = V_child[0,:,:]
qf=ax.quiver(lon_child[::skip,::skip], lat_child[::skip,::skip], U_child[0,::skip,::skip], V_child[0,::skip,::skip], scale=8, transform=ccrs.PlateCarree())
fig.tight_layout()
qk=ax.quiverkey(qf,0.7, 0.7, -0.5, '0.5 m s$^{-1}$', coordinates = 'figure', transform=ccrs.PlateCarree())
# show the time stamp
time_lon=18.6; time_lat=-34.00
tx=ax.text(time_lon, time_lat, str(dates[0]),
horizontalalignment='left',
transform=ccrs.PlateCarree())
# animation function
def animate(i):
#cont.remove()
Uchd = U_child[i,::skip,::skip]
Vchd = V_child[i,::skip,::skip]
qf.set_UVC(Uchd, Vchd)
tx.set_text(str(dates[i]))
anim = FuncAnimation(
fig, animate, interval=200, frames=range(0,240,1)) #len(time)-48
anim.save('Forecast_UV_child_test_processed.gif', writer='imagemagick')
# %% Temperature animations
levels = np.linspace(10,22,20,endpoint=True)
skip = 2
# # setup figure object
fig = plt.figure(figsize=(12,8))
#ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
ax = plt.axes(projection=ccrs.PlateCarree())
#plt.contourf(lon_child, lat_child, SST_child[0,:,:])
ax.set_extent([np.nanmin(lon_child),np.nanmax(lon_child),np.nanmin(lat_child),np.max(lat_child)])
ax.add_feature(cfeature.LAND, zorder=0, edgecolor='black')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle=':')
gl.left_labels = False
gl.top_labels = False
# add the vectors
SST = SST_child[0,:,:]
cont = ax.contourf(lon_child, lat_child, SST,cmap='RdYlBu_r');
cont1 = ax.contour(lon_child, lat_child, SST,levels,colors='k');
ax.contourf(lon_child, lat_child, SST_child[0,:,:],cmap='RdYlBu_r')
cont3 = ax.contour(lon_child, lat_child, SST_child[0,:,:],levels,colors='k')
fig.tight_layout()
# show the time stamp
time_lon=18.6; time_lat=-34.00
tx=ax.text(time_lon, time_lat, str(dates[0]),
horizontalalignment='left',
transform=ccrs.PlateCarree())
# animation function
def animate(i):
global cont
for c in cont.collections:
c.remove() # removes only the contours, leaves the rest intact
#cont.remove()
SST = SST_child[i,:,:]
cont = ax.contourf(lon_child, lat_child, SST,levels,cmap='RdYlBu_r');
global cont1
for c in cont1.collections:
c.remove() # removes only the contours, leaves the rest intact
#cont.remove()
cont1 = ax.contour(lon_child, lat_child, SST,levels,colors='k')
global cont3
for c in cont3.collections:
c.remove() # removes only the contours, leaves the rest intact
#cont.remove()
cont3 = ax.contour(lon_child, lat_child, SST_child[i,:,:],levels,colors='k')
ax.contourf(lon_child, lat_child, SST_child[i,:,:],levels,cmap='RdYlBu_r')
tx.set_text(str(dates[i]))
return cont
anim = FuncAnimation(
fig, animate, interval=200, frames=range(0,240,1)) #len(time)-48
anim.save('Forecast_test_processed.gif', writer='imagemagick')
It would be nice to be able to download python (and other language) generated assets using the CROCO output. For example giffs. To do this the website needs to implement some UI controls for requesting these assets (and generating them), as well as a mechanism for running scripts that may take too long for the duration of a single request (email notifications, or simple file management for logged in users, etc)