{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy\n",
"import os\n",
"import datetime\n",
"import scipy\n",
"from numpy.random import randn\n",
"import numpy as np\n",
"\n",
"from IPython.core.display import display, HTML\n",
"display(HTML(\"\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load and select data from SEAScope\n",
"**Select the surface velocity from GlobCurrent you want to play with (current norm raster)**
\n",
"**Draw a polygon on SEASCope**
\n",
"**Click on the extract button on the polygon window**
\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load data directly from viewer memory\n",
"from SEAScope.lib import get_extracted_data\n",
"extractions = get_extracted_data()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#from SEAScope.lib import load_pyo\n",
"#data_path = ''\n",
"#extractions = load_pyo(data_path)\n",
"## Print a list of path of extracted granules \n",
"print('\\n'.join(extractions.keys()))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Print selected granules**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Play with the data\n",
"** Compute trajectory of fictive swell ray at initial position (x0, y0) using the dispertion relation in presence of surface current :
\n",
"$ \\omega_c = \\omega+ \\vec{k}.\\vec{u}$
\n",
"$ \\frac{dx}{dt}=\\frac{d\\omega_c}{dk}$
\n",
"$ \\frac{dk}{dt}=-\\frac{d\\omega_c}{dx}$
**\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** Fill variable name for your selected velocity **"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Get longitudes, latitudes from gcps and time from first granule\n",
"granule_uri = next( v for i, v in enumerate(extractions.keys()) if i == 0)\n",
"extraction = extractions[granule_uri]\n",
"start = extraction['meta']['start']\n",
"print(extraction['meta']['fields'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## adapt name_var_u and name_var_v to the extracted field names"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"name_var_u = 'ugos'\n",
"name_var_v = 'vgos'\n",
"# name_var_u = 'u'\n",
"# name_var_v = 'v'\n",
"\n",
"extraction = extractions[granule_uri] \n",
"u = extraction['data'][name_var_u]\n",
"v = extraction['data'][name_var_v]\n",
"\n",
"from SEAScope.lib.utils import get_lonlat\n",
"lon2D, lat2D = get_lonlat(extraction, numpy.shape(u))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"velocity = (u**2+v**2)**(1./2)\n",
"print(velocity.shape)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"fig = plt.figure(figsize = (16,10))\n",
"plt.imshow(numpy.flipud(velocity),interpolation='bicubic',cmap='jet', vmin=0 , vmax=2)\n",
"cbar=plt.colorbar()\n",
"X, Y = numpy.meshgrid(numpy.arange(0, u.shape[1]), numpy.arange(0, u.shape[0]))\n",
"Q = plt.quiver(X[::15, ::15], Y[::15, ::15],numpy.clip(numpy.flipud(u[::15,::15]),-2,2), numpy.clip(numpy.flipud(v[::15,::15]),-2,2))\n",
"qk = plt.quiverkey(Q, 0.5, 0.98, 2, r'$2 \\frac{m}{s}$', labelpos='W',\n",
" fontproperties={'weight': 'bold'})\n",
"#plt.xlim(40, 150)\n",
"#plt.ylim(180, 120)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = (lat2D[1,0] - lat2D[0,0])*100000\n",
"res = (lat2D[int(u.shape[0]/2+1), int(u.shape[1]/2)]\n",
" - lat2D[int(u.shape[0]/2), int(u.shape[1]/2)]) * 100000\n",
"print(res)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" ** Play with nstep (number of time steps) and nseed (number of rays) ** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tmpu = extraction['data'][name_var_u]\n",
"urand = numpy.multiply(tmpu,0.02*randn(u.shape[0],u.shape[1]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"def ev(var, x, y):\n",
" if numpy.isnan(x) or numpy.isnan(y):\n",
" return numpy.nan\n",
" nx, ny = numpy.shape(var)\n",
" i = int(x)\n",
" j = int(y)\n",
" if (i > nx - 2) or (j > ny - 2):\n",
" return numpy.nan\n",
" dx = x - i\n",
" dy = y - j\n",
" uix = dx * var[i + 1, j] + (1 - dx) * var[i, j]\n",
" uiy = dx * var[i + 1, j + 1] + (1 - dx) * var[i, j + 1]\n",
" ui = dy * uiy + (1 - dy) * uix\n",
" return ui\n",
"\n",
"\n",
"# Initialize wave ray position\n",
"\n",
"granule_id = None\n",
"# Time step for advection (in sec)\n",
"s=1\n",
"dt = 3600/3/s\n",
"# Number of seeds\n",
"nseed = 120*2 #! Change this\n",
"# Wave number\n",
"\n",
"#10 aout 2019\n",
"wl = 1.56*11**2\n",
"tetad = 25\n",
"c0 = 10 # Starting column\n",
"# Starting line\n",
"l0 = u.shape[0]-150\n",
"# Total time of advection \n",
"nstep = int(150*s) #! Change this\n",
"\n",
"# 12 Aout 2018 16:30\n",
"wl = 350 \n",
"tetad = 50\n",
"c0 = 10 # Starting column\n",
"# Starting line\n",
"l0 = u.shape[0]-150\n",
"# Total time of advection \n",
"nstep = int(100*s) #! Change this\n",
"\n",
"k = 2*numpy.pi/wl\n",
"# Wave direction in radians\n",
"teta0 = tetad*numpy.pi/180\n",
"# Starting column\n",
"\n",
"u = extraction['data'][name_var_u]\n",
"v = extraction['data'][name_var_v]\n",
"res = (lat2D[int(u.shape[0]/2+1), int(u.shape[1]/2)]\n",
" -lat2D[int(u.shape[0]/2), int(u.shape[1]/2)]) * 100000\n",
"x = numpy.full((nstep+1, nseed), numpy.nan)\n",
"y = numpy.full((nstep+1, nseed), numpy.nan)\n",
"lon = numpy.full((nstep+1, nseed), numpy.nan)\n",
"lat = numpy.full((nstep+1, nseed), numpy.nan)\n",
"for c in range(nseed):\n",
" x[0,c] = -l0 - 1.5*c/3\n",
" y[0,c] = c0 + 6*c/3\n",
" lon[0,c] = ev(lon2D, x[0,c], y[0,c])\n",
" lat[0,c] = ev(lat2D, x[0,c], y[0,c])\n",
" teta = teta0\n",
" kx = k*numpy.sin(teta)\n",
" ky = k*numpy.cos(teta)\n",
" for i in range(nstep):\n",
" if x[i, c] is numpy.ma.masked or y[i, c] is numpy.ma.masked:\n",
" continue\n",
"# if (x[i, c]>(u.shape[0]-1)) or (y[i, c]>(u.shape[1]-1)):\n",
"# continue\n",
"# if (x[i, c]<0) or (y[i, c]<0):\n",
"# continue\n",
" omega = numpy.sqrt(9.81*(kx**2+ky**2)) + (ky*ev(u,x[i,c],y[i,c])\n",
" +kx*ev(v,x[i,c],y[i,c]))\n",
" dx = dt*(ev(v,x[i,c],y[i,c])+numpy.sqrt(9.81/(numpy.sqrt(kx**2+ky**2)))\n",
" *numpy.sin(teta)/2)/res\n",
" dy = dt*(ev(u,x[i,c],y[i,c])+numpy.sqrt(9.81/(numpy.sqrt(kx**2+ky**2)))\n",
" *numpy.cos(teta))/2/res\n",
" teta = numpy.arctan2(kx,ky)\n",
" if dx is numpy.ma.masked or dy is numpy.ma.masked:\n",
" continue\n",
" if (x[i, c]+dx<(1-u.shape[0])) or (y[i, c]+dy>(u.shape[1]-1)):\n",
" continue\n",
" if (x[i, c]+dx>0) or (y[i, c]+dy<0):\n",
" continue\n",
" domegax = omega - (numpy.sqrt(9.81*(kx**2+ky**2))\n",
" + (ky*ev(u,x[i,c]+dx,y[i,c])\n",
" + kx*ev(v,x[i,c]+dx,y[i,c])))\n",
" domegay = omega - (numpy.sqrt(9.81*(kx**2+ky**2))\n",
" + (ky*ev(u,x[i,c],y[i,c]+dy)\n",
" + kx*ev(v,x[i,c],y[i,c]+dy)))\n",
"\n",
" x[i+1,c] = x[i,c] + dx\n",
" y[i+1,c] = y[i,c] + dy\n",
" lon[i+1,c] = ev(lon2D, x[i+1,c], y[i+1,c])\n",
" lat[i+1,c] = ev(lat2D, x[i+1,c], y[i+1,c])\n",
" dkx = dt*domegax/dx/res\n",
" dky = dt*domegay/dy/res\n",
" kx = kx + dkx\n",
" ky = ky + dky"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = plt.figure(figsize = (18,10))\n",
"plt.imshow(numpy.flipud(velocity),interpolation='bicubic',cmap='jet', vmin=0 , vmax=2)\n",
"cbar=plt.colorbar()\n",
"Q = plt.quiver(X[::15, ::15], Y[::15, ::15],\n",
" numpy.clip(numpy.flipud(u[::15,::15]),-2,2),\n",
" numpy.clip(numpy.flipud(v[::15,::15]),-2,2))\n",
"qk = plt.quiverkey(Q, 0.5, 0.98, 2, r'$2 \\frac{m}{s}$', labelpos='W',\n",
" fontproperties={'weight': 'bold'})\n",
"for c in range(nseed):\n",
" plt.plot(y[:,c],-x[:,c],color='w',linewidth=1.5)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Export wave rays to SEAScope"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Define Collection:**
\n",
"** A granule has to belong to a collection **
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from SEAScope.lib.utils import create_collection, init_ids\n",
"# Make sure to keep the collection identifier, we will need it later\n",
"collection_id, collection = create_collection('Simulated wave rays')\n",
"print(collection_id)\n",
"\n",
"# IP address and port used to reach the SEAScope standalone application\n",
"host = '127.0.0.1'\n",
"port = 11155\n",
"\n",
"import SEAScope.upload\n",
"with SEAScope.upload.connect(host, port) as link:\n",
"# Send create collection job\n",
" SEAScope.upload.collection(link, collection)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Upload granule to SEAScope:
\n",
"Define shell of granule, create granule
\n",
"The granule will contain each matrix of modified data as what we call a \"field\"
**\n",
"** Enter a name for each field (field_name_time) **\n",
"**Upload collection and granule to SEAScope**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from SEAScope.lib.utils import create_granule\n",
"# Shell of Granule\n",
"# Each position of the trajectory will be used as GCP\n",
"with SEAScope.upload.connect(host, port) as link: \n",
" for c in range(nseed):\n",
" gcps = []\n",
" count = 0\n",
" for i in range(0, nstep+1):\n",
" gcp = {'lon': lon[i,c], 'lat': lat[i,c], 'i': i, 'j': 0}\n",
" if lat[i,c] != 0:\n",
" gcps.append(gcp)\n",
" count=count+1\n",
"\n",
" # Create granule\n",
" stop_dt = start + datetime.timedelta(days=2)\n",
" granule_id, granule = create_granule(collection_id, gcps, start, stop_dt)\n",
"\n",
" # Set variables\n",
" # Enter a name for each variable field\n",
" from SEAScope.lib.utils import set_field\n",
" field_name_time ='time rays'\n",
" set_field(granule, field_name_time, numpy.arange(count)*dt/3600)\n",
" \n",
" # Send create granule job\n",
" SEAScope.upload.granule(link, granule)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Telling SEAScope what to display:
\n",
" The next step is variables definition. **
\n",
"** Enter a name for each variable you want to display (variable_name_time, ...) **"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from SEAScope.lib.utils import create_variable\n",
"\n",
"# Create a variable for each parameter\n",
"variable_name_time = 'time_rays'\n",
"var1 = create_variable(collection, variable_name_time, [field_name_time], dims=1)\n",
"\n",
"# Upload to SEAScope\n",
"with SEAScope.upload.connect(host, port) as link:\n",
" # Send create variable jobs\n",
" SEAScope.upload.variable(link, var1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** Telling SEAScope how to display variables:
\n",
"The rendering options are defined in the following cell.
\n",
" You can play with min, max values, the colormap ... and check how the rendering changes on SEAScope \n",
"**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Rendering configuration for the time\n",
"rcfg1 = var1['rendering']\n",
"rcfg1['min'] = 1\n",
"rcfg1['max'] = 40\n",
"rcfg1['colormap'] = ''\n",
"rcfg1['color'] = [255,255,255]\n",
"rcfg1['opacity'] = 0.25\n",
"rcfg1['zindex'] = .923\n",
"\n",
"\n",
"# Upload to SEAScope\n",
"with SEAScope.upload.connect(host, port) as link:\n",
" # Send create rendering config jobs\n",
" SEAScope.upload.rendering_config(link, rcfg1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** Go back to SEAScope and visually compare your wave rays trajectory with Sentinel2 RGB
\n",
"play with range, transparency ...**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}