{ "cells": [ { "cell_type": "markdown", "id": "99703d03-5db7-4850-add9-77033d4c6217", "metadata": {}, "source": [ "# Generic PCM Photochemistry postprocessing and visualization demonstrator" ] }, { "cell_type": "markdown", "id": "a6b7b35c-fb19-4bda-81dd-ee03df1e4ef8", "metadata": {}, "source": [ "This Notebook will show you how to use the Generic PCM photochemistry postprocessing library and how to make interactive visualization with it. For it to work, you'll need to copy the *photochemistry_postprocessing.py* file along this notebook in the directory containing the output file (*diagfi.nc*) as well as the reaction network file chemnetwork/reactfile (to become *reaction.def*)." ] }, { "cell_type": "markdown", "id": "ed55c2f3-5fa8-481c-b9d8-8e31f93f7993", "metadata": {}, "source": [ "## Loading simulation and calculating reaction rates" ] }, { "cell_type": "markdown", "id": "50a0abe3-e224-4b07-95dc-a424dd35cdcc", "metadata": {}, "source": [ "If you do not have a simulation with photochemistry at hand, you can download an example file with: \n", "```\n", "wget https://web.lmd.jussieu.fr/~mmaurice/photochem_example.zip\n", "unzip photochem_example.zip\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "cd28bac5-65f6-464b-9c2b-6cba1bde9472", "metadata": {}, "outputs": [], "source": [ "import photochem_postproc as pcpp\n", "\n", "sim_path = './photochem_example'\n", "NetCDF_filename = 'diagfi.nc'\n", "\n", "# The simu class is just a wrapper for xr.Dataset\n", "my_sim = pcpp.GPCM_simu(sim_path,NetCDF_filename)" ] }, { "cell_type": "markdown", "id": "703aaead-59d4-4a93-8bfd-25ad8478dee0", "metadata": {}, "source": [ "Now let's load the chemical network. Notice that unlike in this example, the reaction network is normally located in a ***chemnetwork/reactfile*** file." ] }, { "cell_type": "code", "execution_count": null, "id": "355df458-41a8-4d48-a21c-49a427f94ea5", "metadata": {}, "outputs": [], "source": [ "my_sim.network = pcpp.network.from_file(my_sim.path+'/network.def')\n", "my_sim.background_species = 'co2'" ] }, { "cell_type": "markdown", "id": "63d2fe39-aedc-45b9-bfc8-d7670336a876", "metadata": {}, "source": [ "Some reactions' rates are hard-coded and need to be added manually (you should find their rates in *reaction_rate_lib.py*). To do that we first need to define a new reaction and call again **compute_rates** with the new reaction as second argument:" ] }, { "cell_type": "code", "execution_count": null, "id": "bf9f8c56-6615-4f9f-acae-c788a9543770", "metadata": {}, "outputs": [], "source": [ "# First load the parametrization for its rate\n", "from reaction_rate_lib import k_CO_OH_to_CO2_H_JPL_2015\n", "\n", "# Then create the reaction objet (here for the reaction co + oh -> co2 + h):\n", "hard_coded_reaction = pcpp.reaction(['co','oh'],['co2','h'],k_CO_OH_to_CO2_H_JPL_2015)\n", "\n", "my_sim.network.append(hard_coded_reaction)" ] }, { "cell_type": "markdown", "id": "5b193a59-9338-487e-977e-35556a7a405a", "metadata": {}, "source": [ "We are now ready to compute the rates (be careful with memory, as it will add 2 4D fields per chemical species, and 2 4D fields per reaction!)" ] }, { "cell_type": "code", "execution_count": null, "id": "776388e6-0506-428b-a99c-f02283bad4d8", "metadata": {}, "outputs": [], "source": [ "my_sim.compute_rates()" ] }, { "cell_type": "markdown", "id": "042e4416-5ad3-4d1d-bbe7-896619253146", "metadata": {}, "source": [ "## Now let's do some visualization" ] }, { "cell_type": "markdown", "id": "d5e0e556-863a-4f39-b6be-ebe1402a5f95", "metadata": {}, "source": [ "### Static visualization\n", "Here we use the built-in visualization methods of the *simu* class." ] }, { "cell_type": "code", "execution_count": null, "id": "f34cbcfa-339b-46dd-b33e-4c699a118f60", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "779d0619-582e-4628-8161-5c0d9d943261", "metadata": {}, "source": [ "#### Meridional slice" ] }, { "cell_type": "code", "execution_count": null, "id": "f54ea828-ea01-4f8c-9349-e5f051ef7043", "metadata": {}, "outputs": [], "source": [ "plt.subplot(131)\n", "my_sim.plot_meridional_slice('h2o_vap',logcb=True)\n", "plt.title('Time- and longitude-average')\n", "\n", "plt.subplot(132)\n", "my_sim.plot_meridional_slice('h2o_vap',logcb=True,t=0)\n", "plt.title('Time = 0 and longitude-average')\n", "\n", "plt.subplot(133)\n", "my_sim.plot_meridional_slice('h2o_vap',logcb=True,t=0,lon=0)\n", "plt.title('Time = 0 and longitude = 0°')\n", "\n", "plt.subplots_adjust(right=2)" ] }, { "cell_type": "markdown", "id": "3d691c4b-a690-4644-a00a-1ee77c85658a", "metadata": {}, "source": [ "#### Time evolution" ] }, { "cell_type": "code", "execution_count": null, "id": "2328745a-55e1-40d9-86c7-41f04bea872c", "metadata": {}, "outputs": [], "source": [ "plt.subplot(131)\n", "my_sim.plot_time_evolution('h2o_vap',logcb=True)\n", "plt.title('latitude- and longitude-average')\n", "\n", "plt.subplot(132)\n", "my_sim.plot_time_evolution('h2o_vap',logcb=True,lat=0)\n", "plt.title('Equatorial value, longitude-average')\n", "\n", "plt.subplot(133)\n", "my_sim.plot_time_evolution('h2o_vap',logcb=True,lat=0,lon=0)\n", "plt.title('Equatorial value longitude = 0°')\n", "\n", "plt.subplots_adjust(right=2)" ] }, { "cell_type": "markdown", "id": "f287c71d-83eb-4bad-a11e-dfb34e679adf", "metadata": {}, "source": [ "#### Vertical profiles" ] }, { "cell_type": "code", "execution_count": null, "id": "4f544843-012d-41a6-8d7a-c4b40f45a5e1", "metadata": {}, "outputs": [], "source": [ "my_sim.plot_profile('h2o_vap',logx=True,label='Time- and grid-averaged')\n", "my_sim.plot_profile('h2o_vap',t=0,logx=True,label='Time=0, grid-averaged')\n", "my_sim.plot_profile('h2o_vap',t=0,lat=0,logx=True,label='Time=0, equatorial value zonally-averaged')\n", "my_sim.plot_profile('h2o_vap',t=0,lat=0,lon=0,logx=True,label='Time=0, equatorial value, lon = 0°')\n", "plt.legend()\n", "plt.grid()" ] }, { "cell_type": "markdown", "id": "e73804b8-98f3-4283-8c5d-153465db88e7", "metadata": {}, "source": [ "## Interactive visualization\n", "Now we go fancy!" ] }, { "cell_type": "markdown", "id": "7f745af6-fe37-4075-b433-abb03ad3492e", "metadata": {}, "source": [ "#### Define widgets\n", "Let's define some widgets (see [here](https://ipywidgets.readthedocs.io/en/latest/) for documenttion on jupyter's widgets):" ] }, { "cell_type": "code", "execution_count": null, "id": "56b07967-3900-4454-ac0d-29cfae7cc1f9", "metadata": {}, "outputs": [], "source": [ "import ipywidgets as widgets\n", "\n", "# Coordinates\n", "w_lat = widgets.FloatSlider(min=-90, max=90, step=1, description=\"latitude\")\n", "w_lon = widgets.FloatSlider(min=-180, max=180, step=1, description=\"longitude\")\n", "w_alt = widgets.FloatSlider(min=0, max=max(my_sim.data[\"altitude\"]), step=1, description=\"altitude\")\n", "w_time = widgets.FloatSlider(min=0, max=max(my_sim[\"Time\"]), step=0.1, description=\"time\")\n", "\n", "# Fields\n", "w_single_sp = widgets.Select(options=my_sim.network.species, value=\"h2o_vap\", description=\"species\")\n", "w_multiple_sp = widgets.SelectMultiple(options=my_sim.network.species, value=[\"h2o_vap\"], description=\"species\")\n", "w_reactions = widgets.SelectMultiple(options=my_sim.network.reactions.keys(), value=[\"co2 + hv -> co + o\"], description=\"reactions\")\n", "\n", "# Miscelaneous\n", "w_average = widgets.Checkbox(description='show average')" ] }, { "cell_type": "markdown", "id": "ddeec0ba-12c2-4bd9-a38f-98a0e1949f3b", "metadata": {}, "source": [ "#### OH meridional slice at various longitudes\n", "OH is a photolysis product with a short lifetime, so it exist only on the dayside. Scrolling through the longitudes will exhibit this dichotomy." ] }, { "cell_type": "code", "execution_count": null, "id": "56e81d82-fd55-464c-a746-66cf23822957", "metadata": {}, "outputs": [], "source": [ "# Slice plotting unction\n", "def make_slice(lon):\n", " my_sim.plot_meridional_slice('oh',t=0,lon=lon,logcb=True)\n", "\n", "# Define interactive output\n", "out = widgets.interactive_output(make_slice,{'lon':w_lon})\n", "\n", "# Build the output frame\n", "widgets.VBox([w_lon,out])" ] }, { "cell_type": "markdown", "id": "3379c31a-95c6-4ab3-b178-20b706eed1f3", "metadata": {}, "source": [ "#### Water vapor atlas\n", "We can have several action widgets" ] }, { "cell_type": "code", "execution_count": null, "id": "fd7b4103-0436-4bda-bb39-96666c39f332", "metadata": {}, "outputs": [], "source": [ "# Atlas plotting function\n", "def make_atlas(t,alt):\n", " my_sim.plot_atlas('h2o_vap',t=t,alt=alt)\n", " plt.title('t='+str(int(t))+' sol, altitude='+str(int(alt))+' km')\n", "\n", "# Define interactive output\n", "out = widgets.interactive_output(make_atlas,{'t':w_time,'alt':w_alt})\n", "\n", "# Build the output frame\n", "widgets.VBox([w_time,w_alt,out])" ] }, { "cell_type": "markdown", "id": "4b6bb2f8-1aff-4872-9580-7c4bd9598c31", "metadata": {}, "source": [ "#### Temperature profile at various times and locations" ] }, { "cell_type": "code", "execution_count": null, "id": "e4691cae-637b-4555-ac87-d556521a4c3f", "metadata": {}, "outputs": [], "source": [ "# Profile plotting unction\n", "def make_prof(t,lon,lat,avg):\n", " my_sim.plot_profile('temp',t=t,lon=lon,lat=lat,label='lon='+str(int(lon))+'°, lat='+str(int(lat))+'°')\n", " if avg:\n", " my_sim.plot_profile('temp',t=t,label='average')\n", " plt.legend()\n", " plt.grid()\n", "\n", "# Define interactive output\n", "out = widgets.interactive_output(make_prof,{'t':w_time,'lon':w_lon,'lat':w_lat,'avg':w_average})\n", "\n", "# Build the output frame\n", "widgets.HBox([widgets.VBox([w_time,w_lon,w_lat,w_average]),out])" ] }, { "cell_type": "markdown", "id": "9ad5a23a-7ed6-4f56-abbf-3cf16daa087a", "metadata": {}, "source": [ "#### Extensive species visualizer\n", "Combining the above examples for arbitrary species" ] }, { "cell_type": "code", "execution_count": null, "id": "e4db2d8b-6183-4fbe-8ca1-940ef15aaa28", "metadata": {}, "outputs": [], "source": [ "def make_visualizer(sp,t,alt):\n", " \n", " plt.subplot(131) # zonal average\n", " my_sim.plot_meridional_slice(sp,t=t,logcb=True)\n", " plt.plot([-90,90],[alt,alt],ls='--',lw=3,c='white')\n", "\n", " plt.subplot(132) # atlas\n", " my_sim.plot_atlas(sp,t=t,alt=alt)\n", " plt.title('t='+str(int(t))+' sol, altitude='+str(int(alt))+' km')\n", "\n", " plt.subplot(133) # profile\n", " my_sim.plot_profile('temp',t=t)\n", " plt.grid()\n", "\n", " plt.subplots_adjust(right=2)\n", "\n", "out = widgets.interactive_output(make_visualizer,{'sp':w_single_sp,'t':w_time,'alt':w_alt})\n", "\n", "widgets.HBox([widgets.VBox([w_single_sp,w_time,w_alt]),out])" ] }, { "cell_type": "markdown", "id": "b8d39b59-07f2-4988-8e5c-558221c825a1", "metadata": {}, "source": [ "#### Multi-species profiles\n", "Shift+click to select multiple species (Command+click on Mac)" ] }, { "cell_type": "code", "execution_count": null, "id": "9dd59a1c-58c4-41f0-9730-9e97d2607c6a", "metadata": {}, "outputs": [], "source": [ "cmap = plt.get_cmap(\"tab10\")\n", "def make_sp_prof(sps,t,lon,lat,avg):\n", "\n", " for i,sp in enumerate(sps):\n", " my_sim.plot_profile(sp,t=t,lon=lon,lat=lat,logx=True,label=sp,c=cmap(i))\n", " if avg:\n", " my_sim.plot_profile(sp,t=t,logx=True,c=cmap(i),ls='--')\n", " if avg: # just for the legend\n", " plt.plot([],[],c='k',label='lon='+str(int(lon))+'°, lat='+str(int(lat))+'°')\n", " plt.plot([],[],ls='--',c='k',label='average')\n", " \n", " plt.legend()\n", " plt.grid()\n", "\n", "out = widgets.interactive_output(make_sp_prof,{'sps':w_multiple_sp,'t':w_time,'lon':w_lon,'lat':w_lat,'avg':w_average})\n", "\n", "widgets.HBox([widgets.VBox([w_multiple_sp,w_time,w_lon,w_lat,w_average]),out])" ] }, { "cell_type": "markdown", "id": "5fb2c987-345d-4a74-b4bf-fb368eeaac2b", "metadata": {}, "source": [ "#### Species-specific reaction rates" ] }, { "cell_type": "code", "execution_count": null, "id": "ff21a3dc-d44a-4f0e-9aa4-211741bb592d", "metadata": {}, "outputs": [], "source": [ "def make_reaction_rate_viz(sp,t):\n", "\n", " for r in my_sim.network.get_subnetwork({'species':[sp]}):\n", " if sp in r.products:\n", " my_sim.plot_profile('rate ('+r.formula+')',t=t,logx=True,label=r.formula)\n", " elif sp in r.reactants:\n", " my_sim.plot_profile('rate ('+r.formula+')',t=t,logx=True,ls='--',label=r.formula)\n", "\n", " plt.legend()\n", " plt.grid()\n", "\n", "out=widgets.interactive_output(make_reaction_rate_viz,{'sp':w_single_sp,'t':w_time})\n", "\n", "widgets.HBox([widgets.VBox([w_single_sp,w_time]),out])" ] }, { "cell_type": "markdown", "id": "7e0b725d-58f5-4b4e-b170-125b4707df03", "metadata": {}, "source": [ "#### Profile with atlas locator\n", "Here the atlas shows the column mass" ] }, { "cell_type": "code", "execution_count": null, "id": "b5b73171-0101-4656-b2ce-7e398070ebbb", "metadata": {}, "outputs": [], "source": [ "def make_sp_prof_atlas(sps,t,lon,lat,avg):\n", "\n", " plt.subplot(121) # Vertical profile\n", " for i,sp in enumerate(sps):\n", " my_sim.plot_profile(sp,t=t,lon=lon,lat=lat,logx=True,label=sp,c=cmap(i))\n", " if avg:\n", " my_sim.plot_profile(sp,t=t,logx=True,c=cmap(i),ls='--')\n", " if avg: # just for the legend\n", " plt.plot([],[],c='k',label='lon='+str(int(lon))+'°, lat='+str(int(lat))+'°')\n", " plt.plot([],[],ls='--',c='k',label='average')\n", "\n", " plt.legend()\n", " plt.grid()\n", "\n", " plt.subplot(122) # Atlas\n", " my_sim.plot_atlas(sp+'_col',t=t)\n", " plt.scatter([lon],[lat],marker='o',s=[100],c=['tab:red'])\n", "\n", " plt.subplots_adjust(right=2)\n", "\n", "out = widgets.interactive_output(make_sp_prof_atlas,{'sps':w_multiple_sp,'t':w_time,'lon':w_lon,'lat':w_lat,'avg':w_average})\n", "\n", "widgets.HBox([widgets.VBox([w_multiple_sp,w_time,w_lon,w_lat,w_average]),out])" ] }, { "cell_type": "markdown", "id": "239adb94-d3b7-4a9a-81fc-297e72e63639", "metadata": {}, "source": [ "#### Species-specific reaction rates with atlas locator" ] }, { "cell_type": "code", "execution_count": null, "id": "c849fda1-3969-4709-bb17-fdcaff5cbf86", "metadata": {}, "outputs": [], "source": [ "def make_sp_rate_atlas(sp,t,lon,lat,avg):\n", " \n", " plt.subplot(121) # Vertical profile\n", " for r in my_sim.network:\n", " if sp in r.products:\n", " my_sim.plot_profile('rate ('+r.formula+')',t=t,lon=lon if not avg else 'avg',lat=lat if not avg else 'avg',logx=True,label=r.formula)\n", " elif sp in r.reactants:\n", " my_sim.plot_profile('rate ('+r.formula+')',t=t,lon=lon if not avg else 'avg',lat=lat if not avg else 'avg',logx=True,ls='--',label=r.formula)\n", " \n", " plt.legend()\n", " plt.grid()\n", "\n", " plt.subplot(122) # Atlas\n", " my_sim.plot_atlas(sp+'_col',t=t)\n", " plt.scatter([lon],[lat],marker='o',s=[100],c=['tab:red'])\n", "\n", " plt.subplots_adjust(right=1.5)\n", "\n", "out = widgets.interactive_output(make_sp_rate_atlas,{'sp':w_single_sp,'t':w_time,'lon':w_lon,'lat':w_lat,'avg':w_average})\n", "\n", "widgets.HBox([widgets.VBox([w_single_sp,w_time,w_lon,w_lat,w_average]),out])" ] }, { "cell_type": "markdown", "id": "ff6ad79c-d07b-4d81-83d8-fc02e24ba6d4", "metadata": {}, "source": [ "## Chemical network processing" ] }, { "cell_type": "markdown", "id": "da719183-d81b-4dfd-948a-a275e41ae32a", "metadata": {}, "source": [ "### Working with other formats" ] }, { "cell_type": "markdown", "id": "fd6c578a-7a59-4f1f-b180-bc514d4482a4", "metadata": {}, "source": [ "We can import and export chemical network from and into various formats, not only Generic-PCM style." ] }, { "cell_type": "code", "execution_count": null, "id": "e25768d8-914b-4409-947d-b3aa3f3d9c87", "metadata": {}, "outputs": [], "source": [ "# Let's export our simulation's network into a file readable by the VULCAN model\n", "my_sim.network.to_file('photochem_example/VULCAN_network.txt',format='vulcan')" ] }, { "cell_type": "markdown", "id": "e12f6270-65ba-44d3-8ec7-28292fb50cce", "metadata": {}, "source": [ "Some notifications have been issued because one has to manually add some extra informations for some reactions. But you the file ***photochem_example/VULCAN_network.txt*** should have been correctly written, and you can use it to run a VULCAN simulation once you have added the required informations." ] }, { "cell_type": "markdown", "id": "0b258f20-f574-4382-b394-afcf171d1993", "metadata": {}, "source": [ "### Extracting subnetwork" ] }, { "cell_type": "markdown", "id": "819b6abf-0ec5-4b53-8af5-31da362f6b86", "metadata": {}, "source": [ "We can also use these tools to design and export subnetworks:" ] }, { "cell_type": "code", "execution_count": null, "id": "1e1585ff-dc78-44ae-9e38-3025fde45a59", "metadata": {}, "outputs": [], "source": [ "# Select by type:\n", "print('Photolysis sub-network:')\n", "print('=======================')\n", "for r in my_sim.network.get_subnetwork({'type':[pcpp.photolysis]}):\n", " print(r.formula)\n", "\n", "# Select by species:\n", "print('')\n", "print('HOx sub-network:')\n", "print('=================')\n", "for r in my_sim.network.get_subnetwork({'species':['h','oh','ho2']}):\n", " print(r.formula)\n", "\n", "# Select by element:\n", "print('')\n", "print('C sub-network:')\n", "print('===============')\n", "for r in my_sim.network.get_subnetwork({'elements':['c']}):\n", " print(r.formula)" ] }, { "cell_type": "markdown", "id": "4ea46279-9234-4a53-800d-891efa4ec9d1", "metadata": {}, "source": [ "### Exporting custom subnetwork\n", "You can make your own custom chemical network (ideally selecting reactions from a large database - here we just work with our original simulation's reduced network). You may want to save the associated ***traceur.def*** file (notice that you'll have to add manually parameters, but at least you won't forget tracers.)" ] }, { "cell_type": "code", "execution_count": null, "id": "513f6d4e-a0e2-4095-9a7c-acafc79a465b", "metadata": {}, "outputs": [], "source": [ "from ipywidgets import Layout\n", "\n", "def select_reactions(reactions,net_filename,trac_filename,format):\n", " save_network_button = widgets.Button(description='Save reactfile')\n", " save_tracers_button = widgets.Button(description='Save tracfile')\n", " \n", " def on_save_network_button_clicked(b):\n", " to_save = pcpp.network()\n", " for r in reactions:\n", " to_save.append(my_sim.network[r])\n", " to_save.to_file(net_filename,format=format)\n", " \n", " def on_save_tracers_button_clicked(b):\n", " to_save = pcpp.network()\n", " for r in reactions:\n", " to_save.append(my_sim.network[r])\n", " to_save.save_traceur_file(trac_filename)\n", " \n", " save_network_button.on_click(on_save_network_button_clicked)\n", " save_tracers_button.on_click(on_save_tracers_button_clicked)\n", " display(save_network_button, widgets.Output())\n", " display(save_tracers_button, widgets.Output())\n", " \n", "w_reactions = widgets.SelectMultiple(options=my_sim.network.reactions.keys(),layout=Layout(height='400px'))\n", "w_network_filename = widgets.Text(value='photochem_example/my_custom_network.def',description='network filename:',disabled=False)\n", "w_tracers_filename = widgets.Text(value='photochem_example/my_custom_traceur.def',description='tracers filename:',disabled=False)\n", "w_format = widgets.Select(options=['GPCM','vulcan'],value='GPCM',description='format:')\n", "\n", "out = widgets.interactive_output(select_reactions,{'reactions':w_reactions,'net_filename':w_network_filename,\n", " 'trac_filename':w_tracers_filename,'format':w_format})\n", "widgets.HBox([w_reactions,widgets.VBox([w_format,w_network_filename,w_tracers_filename,out])])" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.7" } }, "nbformat": 4, "nbformat_minor": 5 }