"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Handling 3D Inputs with PICASO\n",
"\n",
"Many GCM groups have their files stored as `netCDF` files. Therefore, you may already be used to `xarray` format. If that is the case you will be able to directly input your xarray formatted data to `PICASO` to get out post-processed spectra. If not though, this tutorial will walk you through how to structure your data in xarray format. \n",
"\n",
"What you will learn: \n",
"\n",
"1. How to convert traditional numpy arrays to `xarray` formatted data, which is common for 3D GCM output\n",
"2. How to regrid using `xarray` and `xesmf`'s `regridder`\n",
"3. How to use `PICASO` built in function (which use #2)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"\n",
"from picaso import justdoit as jdi"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Two new packages you will need, that are not required components of other models: xesmf and xarray. You can read more about the installing these packages here: \n",
"\n",
"[Install XESMF](https://xesmf.readthedocs.io/en/latest/installation.html)\n",
"- will only be needed if you want to use their really handy regridding tools \n",
"\n",
"[Install xarray](http://xarray.pydata.org/en/stable/getting-started-guide/installing.html)\n",
"- needed for all `PICASO` 3d operations"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xesmf as xe\n",
"import xarray as xr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will begin with an example file from the MIT GCM group (courtesy of Tiffany Kataria). "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gcm_out = jdi.HJ_pt_3d()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are going to go through the motions of converting basic `numpy` arrays to `xarray` format. **If you already understand xarrays you may skip to** the [`picaso` section](#regrid-3d-gcm-with-picaso)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gcm_out.keys()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, pressure, temperature, and kzz are all on a coordinate system that is : \n",
"\n",
"n_longitude (128) x n_latitude (64) x n_pressure (53)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gcm_out['temperature'].shape, len(gcm_out['longitude']), len(gcm_out['latitude'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## `xarray` tutorial: Convert `numpy` arrays to `xarray` DataSet\n",
"\n",
"The comments with `required` next to them indicate that they are required for `picaso` to create a spectrum"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# create data\n",
"data = gcm_out['temperature']\n",
"\n",
"# create coords\n",
"lon = gcm_out['longitude']\n",
"lat = gcm_out['latitude']\n",
"pres = gcm_out['pressure'][0,0,:]\n",
"\n",
"# put data into a dataset\n",
"ds = xr.Dataset(\n",
" data_vars=dict(\n",
" temperature=([\"lon\", \"lat\",\"pressure\"], data,{'units': 'Kelvin'})#, required\n",
" #kzz = ([\"x\", \"y\",\"z\"], gcm_out['kzz'])#could add other data components if wanted\n",
" ),\n",
" coords=dict(\n",
" lon=([\"lon\"], lon,{'units': 'degrees'}),#required\n",
" lat=([\"lat\"], lat,{'units': 'degrees'}),#required\n",
" pressure=([\"pressure\"], pres,{'units': 'bar'})#required*\n",
" ),\n",
" attrs=dict(description=\"coords with vectors\"),\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## `xarray` tutorial: Easy plotting"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds['temperature'].isel(pressure=50).plot(x='lon',y='lat')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## `xesmf` tutorial: Step-by-step regrid 3D GCM \n",
"\n",
"The biggest complication with moving to 3D is making sure that the latitude/longitude grids of some users GCM and `picaso` line up properly. `picaso` computes a flux integration on specific `gauss` and `tchebychev` angles. **Meaning, your GCM input will need to be regridded to fit our angles.** Luckily, as you will see below, it is very easy to do this!\n",
"\n",
"First, we will show you how this is done using `xesmf`, then we will introduce the `PICASO` function that leverages these same techniques. \n",
"\n",
"### Step 1) Get latitude/longitude grid used by `picaso`\n",
"\n",
"**G**auss angles are essentially equivalent to lon**G**itudes \n",
"\n",
"**T**chebychev angles are essentially equivalent to la**T**itudes. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n_gauss_angles =10\n",
"n_chebychev_angles=10\n",
"\n",
"gangle, gweight, tangle, tweight = jdi.get_angles_3d(n_gauss_angles, n_chebychev_angles)\n",
"ubar0, ubar1, cos_theta, latitude, longitude = jdi.compute_disco(n_gauss_angles, n_chebychev_angles, gangle, tangle, phase_angle=0)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Step 2) Create the `xesmf` regridder "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds_out = xr.Dataset({'lon': (['lon'], longitude*180/np.pi),\n",
" 'lat': (['lat'], latitude*180/np.pi),\n",
" }\n",
" )\n",
"\n",
"regridder = xe.Regridder(ds, ds_out, 'bilinear')\n",
"ds_out = regridder(ds,keep_attrs=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All done!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds_out['temperature'].isel(pressure=10).plot(x='lon', y ='lat')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Regrid 3D GCM with PICASO\n",
"\n",
"The above code is all the PICASO built in function does -- in addition to completing some checks to make sure that your model run is on the same grid as what you are wanting PICASO run. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds = jdi.HJ_pt_3d(as_xarray=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### DIY Option 1\n",
"\n",
"For completeness, first here is the function `picaso` uses internally. You might use this if you want to manipulate results before supplying `picaso` with your final input."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#regrid yourself\n",
"out_ds = jdi.regrid_xarray(ds, num_gangle=20, \n",
" num_tangle=20, phase_angle=0) \n",
"#then supply to picaso\n",
"case_3d = jdi.inputs()\n",
"case_3d.phase_angle(0, num_gangle=20, num_tangle=20)\n",
"\n",
"#here, regrid is false because you have already done it yourself\n",
"case_3d.atmosphere_3d(out_ds,regrid=False,plot=True,verbose=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Easiest Option 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the regular `PICASO` workflow to regrid your 3d input"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"case_3d = jdi.inputs()\n",
"case_3d.phase_angle(0, num_gangle=20, num_tangle=20)\n",
"#regrid is True as you will do it within the function\n",
"case_3d.atmosphere_3d(ds,regrid=True,plot=True,verbose=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have solved regridding our temperature-pressure profile. However, PICASO has notified us of another step we must take before running a spectrum: *verbose=True;Only one data variable included. Make sure to add in chemical abundances before trying to run spectra.* \n",
"\n",
"In the next notebook you will see how to add chemistry and/or clouds to your 3D input. "
]
}
],
"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.5"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}