Skip to content

Some updates to sailship code #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ammedd merged 2 commits into from
Nov 2, 2023
Merged

Some updates to sailship code #2

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
857ff5e
Some updates to sailship code
erikvansebille Oct 20, 2023
fc6b88a
Using pset_CTD for particlefile creation
erikvansebille Oct 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 47 additions & 40 deletions sailship.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,12 +43,10 @@
"dataset_folder = \"/nethome/0448257/Data\"\n",
"filenames = {\n",
" \"U\": f\"{dataset_folder}/glo12_rg_1d-m_2021010*_3D-uovo_hcst_R20210113.nc\",\n",
" \"Usample\": f\"{dataset_folder}/glo12_rg_1d-m_2021010*_3D-uovo_hcst_R20210113.nc\",\n",
" \"V\": f\"{dataset_folder}/glo12_rg_1d-m_2021010*_3D-uovo_hcst_R20210113.nc\",\n",
" \"Vsample\": f\"{dataset_folder}/glo12_rg_1d-m_2021010*_3D-uovo_hcst_R20210113.nc\",\n",
" \"S\": f\"{dataset_folder}/glo12_rg_1d-m_2021010*_3D-so_hcst_R20210113.nc\",\n",
" \"T\": f\"{dataset_folder}/glo12_rg_1d-m_2021010*_3D-thetao_hcst_R20210113.nc\"}\n",
"variables = {'U': 'uo', 'Usample': 'uo', 'V': 'vo', 'Vsample': 'vo', 'S':'so', 'T':'thetao'}\n",
"variables = {'U': 'uo', 'V': 'vo', 'S':'so', 'T':'thetao'}\n",
"dimensions = {'lon': 'longitude', 'lat': 'latitude', 'time': 'time', 'depth':'depth'}\n",
"\n",
"# create the fieldset and set interpolation methods\n",
Expand All @@ -61,8 +59,9 @@
"bathymetry_variables = ('bathymetry', 'deptho')\n",
"bathymetry_dimensions = {'lon': 'longitude', 'lat': 'latitude'}\n",
"bathymetry_field = Field.from_netcdf(bathymetry_file, bathymetry_variables, bathymetry_dimensions)\n",
"fieldset.add_field(bathymetry_field) \n",
"fieldset.add_field(bathymetry_field)\n",
"fieldset.add_constant('z_start',0.5) #TODO\n",
"fieldset.computeTimeChunk(0,1)\n",
"\n",
"# set initial location #TODO should be input to function or come from leafmap #TODO remove points outside domain\n",
"coords = [[-94.394531, -4.390229], [-106.523437, -0.175781], [-124.541016, -6.053161], [-140.537109, -0.527336]]"
Expand Down Expand Up @@ -97,7 +96,7 @@
" g = pyproj.Geod(ellps='WGS84')\n",
" # r = g.inv_intermediate(startlong, startlat, endlong, endlat, del_s = 1080)\n",
" r = g.inv_intermediate(startlong, startlat, endlong, endlat, initial_idx=0, return_back_azimuth=False, npts=5)\n",
" # store all intermediate points \n",
" # store all intermediate points\n",
" sample_lons.append(r.lons) # stored as a list of arrays\n",
" sample_lats.append(r.lats)\n",
"\n",
Expand All @@ -115,7 +114,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -127,12 +126,12 @@
"\n",
"# define ADCP sampling function from Usample and Vsample (because of A grid)\n",
"def SampleVel(particle, fieldset, time):\n",
" particle.U = fieldset.Usample[time, particle.depth, particle.lat, particle.lon]\n",
" particle.V = fieldset.Vsample[time, particle.depth, particle.lat, particle.lon]\n",
" particle.U, particle.V = fieldset.UV.eval(time, particle.depth, particle.lat, particle.lon, applyConversion=False)\n",
" # particle.V = fieldset.V.eval(time, particle.depth, particle.lat, particle.lon, applyConversion=False)\n",
"\n",
"# Create CTD like particles to sample the ocean\n",
"class CTDParticle(JITParticle):\n",
" \"\"\"Define a new particle class that does CTD like measurements\"\"\" \n",
" \"\"\"Define a new particle class that does CTD like measurements\"\"\"\n",
" salinity = Variable(\"salinity\", initial=np.nan)\n",
" temperature = Variable(\"temperature\", initial=np.nan)\n",
" raising = Variable(\"raising\", dtype=np.int32, initial=0.0)\n",
Expand All @@ -145,23 +144,23 @@
" if particle.raising == 0:\n",
" # Sinking with vertical_speed until near seafloor\n",
" particle_ddepth = vertical_speed * particle.dt\n",
" if particle.depth >= (seafloor - 20): \n",
" if particle.depth >= (seafloor - 20):\n",
" particle.raising = 1\n",
"\n",
" elif particle.raising == 1:\n",
" # Rising with vertical_speed until depth is 1.5 m\n",
" while particle.depth > 1.5:\n",
" particle_ddepth = vertical_speed * particle.dt\n",
" particle_ddepth = -vertical_speed * particle.dt # TODO: should this not be a negative value (EvS)?\n",
" if particle.depth <= 1.5:\n",
" # to break the loop ...\n",
" particle.state = 4\n",
" print(\"CTD cast finished\")\n",
"\n",
"# define function sampling Salinity \n",
"# define function sampling Salinity\n",
"def SampleS(particle, fieldset, time):\n",
" particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon]\n",
"\n",
"# define function sampling Temperature \n",
"# define function sampling Temperature\n",
"def SampleT(particle, fieldset, time):\n",
" particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]"
]
Expand All @@ -175,27 +174,33 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 5,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'ParticleSet' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m/nethome/0448257/Virtual_ship_classroom/sailship.ipynb Cell 10\u001b[0m line \u001b[0;36m4\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=1'>2</a>\u001b[0m depthnum \u001b[39m=\u001b[39m \u001b[39m5\u001b[39m\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a>\u001b[0m \u001b[39m# Initiate ADCP like particle set \u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=3'>4</a>\u001b[0m pset \u001b[39m=\u001b[39m ParticleSet\u001b[39m.\u001b[39mfrom_list(\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a>\u001b[0m fieldset\u001b[39m=\u001b[39mfieldset, pclass\u001b[39m=\u001b[39mADCPParticle, lon\u001b[39m=\u001b[39mnp\u001b[39m.\u001b[39mfull(depthnum,sample_lons[\u001b[39m0\u001b[39m]), lat\u001b[39m=\u001b[39mnp\u001b[39m.\u001b[39mfull(depthnum,sample_lats[\u001b[39m0\u001b[39m]), depth\u001b[39m=\u001b[39mnp\u001b[39m.\u001b[39mlinspace(\u001b[39m5\u001b[39m, \u001b[39m1000\u001b[39m, num\u001b[39m=\u001b[39mdepthnum)\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a>\u001b[0m )\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=7'>8</a>\u001b[0m \u001b[39m# create a ParticleFile to store the ADCP output\u001b[39;00m\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Blorenz-compute/nethome/0448257/Virtual_ship_classroom/sailship.ipynb#X12sdnNjb2RlLXJlbW90ZQ%3D%3D?line=8'>9</a>\u001b[0m adcp_output_file \u001b[39m=\u001b[39m pset\u001b[39m.\u001b[39mParticleFile(name\u001b[39m=\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.results/sailship_adcp.zarr\u001b[39m\u001b[39m\"\u001b[39m)\n",
"\u001b[0;31mNameError\u001b[0m: name 'ParticleSet' is not defined"
"name": "stdout",
"output_type": "stream",
"text": [
"0 P[0](lon=-94.394531, lat=-4.390229, depth=5.000000, U=0.000000, V=0.000000, time=0.000000)\n",
"P[1](lon=-94.394531, lat=-4.390229, depth=253.750000, U=0.000000, V=0.000000, time=0.000000)\n",
"P[2](lon=-94.394531, lat=-4.390229, depth=502.500000, U=0.000000, V=0.000000, time=0.000000)\n",
"P[3](lon=-94.394531, lat=-4.390229, depth=751.250000, U=0.000000, V=0.000000, time=0.000000)\n",
"P[4](lon=-94.394531, lat=-4.390229, depth=1000.000000, U=0.000000, V=0.000000, time=0.000000)\n",
"INFO: Output files are stored in .results/sailship_adcp.zarr.\n",
"100%|██████████| 1.0/1.0 [00:02<00:00, 2.62s/it]\n",
"Writing time 0.0\n",
"CTD time!\n",
"INFO: Output files are stored in CTD_test_1.zarr.\n",
" 23%|██▎ | 4080.0/18000.0 [01:45<04:55, 47.18it/s]"
]
}
],
"source": [
"# Create ADCP like particle set TODO accurate depths \n",
"# Create ADCP like particle set TODO accurate depths\n",
"depthnum = 5\n",
"# Initiate ADCP like particle set \n",
"# Initiate ADCP like particle set\n",
"pset = ParticleSet.from_list(\n",
" fieldset=fieldset, pclass=ADCPParticle, lon=np.full(depthnum,sample_lons[0]), lat=np.full(depthnum,sample_lats[0]), depth=np.linspace(5, 1000, num=depthnum)\n",
" fieldset=fieldset, pclass=ADCPParticle, lon=np.full(depthnum,sample_lons[0]), lat=np.full(depthnum,sample_lats[0]), depth=np.linspace(5, 1000, num=depthnum), time=0\n",
")\n",
"\n",
"# create a ParticleFile to store the ADCP output\n",
Expand All @@ -206,6 +211,7 @@
"\n",
"# initialize CTD station number\n",
"ctd = 0\n",
"ctddt = timedelta(seconds=120) # timestep of CTD output (for speed-up)\n",
"\n",
"# run the model for the length of the sample_lons list\n",
"for i in range(len(sample_lons)):\n",
Expand All @@ -214,35 +220,36 @@
" pset.lon[:] = sample_lons[i]\n",
" pset.lat[:] = sample_lats[i]\n",
"\n",
" # update the particle time\n",
" total_time += timedelta(hours=1).total_seconds()\n",
" pset.time[:] = total_time\n",
" print(i, pset)\n",
"\n",
" # execute the kernels to sample U and V\n",
" pset.execute(SampleVel, dt=0, output_file=adcp_output_file)\n",
" pset.execute(SampleVel, dt=1, runtime=1, output_file=adcp_output_file)\n",
" print(f\"Writing time {pset.time[0]}\")\n",
"\n",
" # check if we are at a CTD station\n",
" if sample_lons[i] == coords[ctd][0] and sample_lats[i] == coords[ctd][1]:\n",
" print('CTD time!')\n",
" ctd += 1\n",
"\n",
" # create a ParticleFile to store the CTD output\n",
" ctd_output_file = pset.ParticleFile(name=f\"CTD_test_{ctd}.zarr\")\n",
"\n",
" # release CTD particle\n",
" pset_CTD = ParticleSet(fieldset=fieldset, pclass=CTDParticle, lon=sample_lons[i], lat=sample_lats[i], depth=fieldset.z_start, time=total_time) \n",
" pset_CTD = ParticleSet(fieldset=fieldset, pclass=CTDParticle, lon=sample_lons[i], lat=sample_lats[i], depth=fieldset.z_start, time=total_time)\n",
"\n",
" # create a ParticleFile to store the CTD output\n",
" ctd_output_file = pset_CTD.ParticleFile(name=f\"CTD_test_{ctd}.zarr\", outputdt=ctddt)\n",
"\n",
" # record the temperature and salinity of the particle\n",
" pset_CTD.execute([CTDcast, SampleS, SampleT], runtime=timedelta(hours=5), dt=timedelta(seconds=1), output_file=ctd_output_file)\n",
" pset_CTD.execute([CTDcast, SampleS, SampleT], runtime=timedelta(hours=5), dt=ctddt, output_file=ctd_output_file)\n",
" total_time = pset_CTD.time[0]\n",
" print(total_time)\n",
"\n"
"\n",
" # update the particle time\n",
" total_time += timedelta(hours=1).total_seconds()\n",
" pset.time[:] = total_time\n"
]
},
{
"cell_type": "code",
"execution_count": 26,
"execution_count": null,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -274,7 +281,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": null,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -340,7 +347,7 @@
},
{
"cell_type": "code",
"execution_count": 28,
"execution_count": null,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -387,7 +394,7 @@
},
{
"cell_type": "code",
"execution_count": 31,
"execution_count": null,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -448,7 +455,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.5"
"version": "3.11.6"
},
"orig_nbformat": 4
},
Expand Down