From 0137a7427dd9bd44746f11fdf03b28ee07d713b0 Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 20:44:53 +0800 Subject: [PATCH 1/8] feat: add `interp1d_vertical_pressure2altitude` --- .github/workflows/docs_readthedocs.yml | 70 +++ .../plot_basic_statistical_analysis.ipynb | 568 ++---------------- .../example/plot_formatting_coordinates.ipynb | 2 +- .../plot_geographic_finite_difference.ipynb | 439 ++------------ docs/example/plot_interp.ipynb | 110 +++- docs/example/plot_taylor_diagram.ipynb | 2 +- docs/example/plot_time_scale_average.ipynb | 2 +- docs/source/dynamic_docs/plot_interp.py | 341 +++++++++++ src/easyclimate/interp/__init__.py | 1 + src/easyclimate/interp/pressure2altitude.py | 97 +++ test/test_interp_pressure2altitude.py | 221 +++++++ 11 files changed, 961 insertions(+), 892 deletions(-) create mode 100644 .github/workflows/docs_readthedocs.yml create mode 100644 src/easyclimate/interp/pressure2altitude.py create mode 100644 test/test_interp_pressure2altitude.py diff --git a/.github/workflows/docs_readthedocs.yml b/.github/workflows/docs_readthedocs.yml new file mode 100644 index 00000000..d10da1df --- /dev/null +++ b/.github/workflows/docs_readthedocs.yml @@ -0,0 +1,70 @@ +name: docs_readthedocs + +# execute this workflow automatically when a we push to master +on: + push: + branches: [ main, dev ] + +# Sets permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages +# https://github.com/ad-m/github-push-action/issues/52 +permissions: write-all + +jobs: + + build_docs_job: + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + - name: Checkout 🛎️ + uses: actions/checkout@v4.2.2 + + - name: Set up Python + uses: actions/setup-python@v5.3.0 + with: + python-version: '3.10' + + - name: Install OptiPNG and Build 🔧 + run: | + sudo apt-get install -y optipng + + - name: Install dependencies 💻 + run: | + pip install -r docs/requirements.txt + + - name: make the sphinx docs + run: | + # make -C docs clean + make -C docs html + + - name: Init new repo in dist folder and commit generated files + run: | + cd docs/build/html/ + git init + touch .nojekyll + git add -A + git config --local user.email "action@github.com" + git config --local user.name "GitHub Action" + git commit -m 'deploy' + + - name: Force push to destination branch 🚀 + uses: ad-m/github-push-action@v0.8.0 + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + branch: gh-pages + force: true + directory: ./docs/build/html + + - name: Upload artifact + uses: actions/upload-pages-artifact@v3 + with: + path: ./docs/build/html + + + - name: Trigger RTDs build + uses: dfm/rtds-action@v1 + with: + webhook_url: ${{ secrets.RTDS_WEBHOOK_URL }} + webhook_token: ${{ secrets.RTDS_WEBHOOK_TOKEN }} + commit_ref: ${{ github.ref }} diff --git a/docs/example/plot_basic_statistical_analysis.ipynb b/docs/example/plot_basic_statistical_analysis.ipynb index d310ccca..bb4d89ce 100644 --- a/docs/example/plot_basic_statistical_analysis.ipynb +++ b/docs/example/plot_basic_statistical_analysis.ipynb @@ -4,10 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "\n", - "# Basic Statistical Analysis\n", - "\n", - "Before proceeding with all the steps, first import some necessary libraries and packages\n" + "\n# Basic Statistical Analysis\n\nBefore proceeding with all the steps, first import some necessary libraries and packages\n" ] }, { @@ -18,20 +15,14 @@ }, "outputs": [], "source": [ - "import easyclimate as ecl\n", - "import matplotlib.pyplot as plt\n", - "import cartopy.crs as ccrs" + "import easyclimate as ecl\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Obtain the sea ice concentration (SIC) data from the Barents-Kara Seas (30°−90°E, 65°−85°N).\n", - "\n", - ".. seealso::\n", - " Luo, B., Luo, D., Ge, Y. et al. Origins of Barents-Kara sea-ice interannual variability modulated by the Atlantic pathway of El Niño–Southern Oscillation. Nat Commun 14, 585 (2023). https://doi.org/10.1038/s41467-023-36136-5\n", - "\n" + "Obtain the sea ice concentration (SIC) data from the Barents-Kara Seas (30\u00b0\u221290\u00b0E, 65\u00b0\u221285\u00b0N).\n\n.. seealso::\n Luo, B., Luo, D., Ge, Y. et al. Origins of Barents-Kara sea-ice interannual variability modulated by the Atlantic pathway of El Ni\u00f1o\u2013Southern Oscillation. Nat Commun 14, 585 (2023). https://doi.org/10.1038/s41467-023-36136-5\n\n" ] }, { @@ -42,19 +33,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea = ecl.open_tutorial_dataset(\"mini_HadISST_ice\").sic\n", - "sic_data_Barents_Sea" + "sic_data_Barents_Sea = ecl.open_tutorial_dataset(\"mini_HadISST_ice\").sic\nsic_data_Barents_Sea" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "And tropical SST dataset.\n", - "\n", - ".. seealso::\n", - " Rayner, N. A.; Parker, D. E.; Horton, E. B.; Folland, C. K.; Alexander, L. V.; Rowell, D. P.; Kent, E. C.; Kaplan, A. (2003) Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century J. Geophys. Res.Vol. 108, No. D14, 4407 https://doi.org/10.1029/2002JD002670 (pdf ~9Mb)\n", - "\n" + "And tropical SST dataset.\n\n.. seealso::\n Rayner, N. A.; Parker, D. E.; Horton, E. B.; Folland, C. K.; Alexander, L. V.; Rowell, D. P.; Kent, E. C.; Kaplan, A. (2003) Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century J. Geophys. Res.Vol. 108, No. D14, 4407 https://doi.org/10.1029/2002JD002670 (pdf ~9Mb)\n\n" ] }, { @@ -65,17 +51,14 @@ }, "outputs": [], "source": [ - "sst_data = ecl.open_tutorial_dataset(\"mini_HadISST_sst\").sst\n", - "sst_data" + "sst_data = ecl.open_tutorial_dataset(\"mini_HadISST_sst\").sst\nsst_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Mean States for Special Month\n", - ":py:func:`easyclimate.get_specific_months_data ` allows us to easily obtain data on the SIC for December alone.\n", - "\n" + "## Mean States for Special Month\n:py:func:`easyclimate.get_specific_months_data ` allows us to easily obtain data on the SIC for December alone.\n\n" ] }, { @@ -86,16 +69,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12 = ecl.get_specific_months_data(sic_data_Barents_Sea, 12)\n", - "sic_data_Barents_Sea_12" + "sic_data_Barents_Sea_12 = ecl.get_specific_months_data(sic_data_Barents_Sea, 12)\nsic_data_Barents_Sea_12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we try to draw the mean states of the SIC in the Barents-Kara for the December.\n", - "\n" + "Now we try to draw the mean states of the SIC in the Barents-Kara for the December.\n\n" ] }, { @@ -106,39 +87,14 @@ }, "outputs": [], "source": [ - "draw_sic_mean_state = sic_data_Barents_Sea_12.mean(dim=\"time\")\n", - "\n", - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "draw_sic_mean_state.plot.contourf(\n", - " ax=ax,\n", - " # projection on data\n", - " transform=ccrs.PlateCarree(),\n", - " # Colorbar is placed at the bottom\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"Blues\",\n", - " levels=21,\n", - ")" + "draw_sic_mean_state = sic_data_Barents_Sea_12.mean(dim=\"time\")\n\nfig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\ndraw_sic_mean_state.plot.contourf(\n ax=ax,\n # projection on data\n transform=ccrs.PlateCarree(),\n # Colorbar is placed at the bottom\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"Blues\",\n levels=21,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Linear Trend\n", - "As we all know, the area of Arctic sea ice has been decreasing more and more in recent years\n", - "due to the impact of global warming. We can obtain the change of SIC by solving\n", - "the linear trend of SIC data from 1981-2022.\n", - ":py:func:`easyclimate.calc_linregress_spatial ` can provide the calculation\n", - "results of solving the linear trend for each grid point.\n", - "\n" + "## Linear Trend\nAs we all know, the area of Arctic sea ice has been decreasing more and more in recent years\ndue to the impact of global warming. We can obtain the change of SIC by solving\nthe linear trend of SIC data from 1981-2022.\n:py:func:`easyclimate.calc_linregress_spatial ` can provide the calculation\nresults of solving the linear trend for each grid point.\n\n" ] }, { @@ -149,18 +105,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_linear_trend = ecl.calc_linregress_spatial(\n", - " sic_data_Barents_Sea_12, dim=\"time\"\n", - ").compute()\n", - "sic_data_Barents_Sea_12_linear_trend" + "sic_data_Barents_Sea_12_linear_trend = ecl.calc_linregress_spatial(\n sic_data_Barents_Sea_12, dim=\"time\"\n).compute()\nsic_data_Barents_Sea_12_linear_trend" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The `slope` is our desired linear trend, let's try to plot the linear trend of each grid point.\n", - "\n" + "The `slope` is our desired linear trend, let's try to plot the linear trend of each grid point.\n\n" ] }, { @@ -171,33 +123,14 @@ }, "outputs": [], "source": [ - "draw_sic_slope = sic_data_Barents_Sea_12_linear_trend.slope\n", - "\n", - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "draw_sic_slope.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")" + "draw_sic_slope = sic_data_Barents_Sea_12_linear_trend.slope\n\nfig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\ndraw_sic_slope.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"RdBu_r\",\n levels=21,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The `pvalue` is the corresponding p-value, and we can determine a significance level (e.g., significance level is set to 0.05)\n", - "in order to plot the region of significance.\n", - "\n" + "The `pvalue` is the corresponding p-value, and we can determine a significance level (e.g., significance level is set to 0.05)\nin order to plot the region of significance.\n\n" ] }, { @@ -208,29 +141,14 @@ }, "outputs": [], "source": [ - "draw_sic_pvalue = sic_data_Barents_Sea_12_linear_trend.pvalue\n", - "\n", - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "ecl.plot.draw_significant_area_contourf(\n", - " draw_sic_pvalue, ax=ax, thresh=0.05, transform=ccrs.PlateCarree()\n", - ")" + "draw_sic_pvalue = sic_data_Barents_Sea_12_linear_trend.pvalue\n\nfig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\necl.plot.draw_significant_area_contourf(\n draw_sic_pvalue, ax=ax, thresh=0.05, transform=ccrs.PlateCarree()\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Further, we can superimpose the linear trend and the region of significance to study the linear trend of\n", - "the region of significance (since the linear trend of the region of non-significance is often spurious).\n", - "\n" + "Further, we can superimpose the linear trend and the region of significance to study the linear trend of\nthe region of significance (since the linear trend of the region of non-significance is often spurious).\n\n" ] }, { @@ -241,51 +159,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "# SIC slope\n", - "draw_sic_slope.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "\n", - "# SIC 95% significant level\n", - "ecl.plot.draw_significant_area_contourf(\n", - " draw_sic_pvalue, ax=ax, thresh=0.05, transform=ccrs.PlateCarree()\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\n# SIC slope\ndraw_sic_slope.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"RdBu_r\",\n levels=21,\n)\n\n# SIC 95% significant level\necl.plot.draw_significant_area_contourf(\n draw_sic_pvalue, ax=ax, thresh=0.05, transform=ccrs.PlateCarree()\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Regression\n", - "Regression analysis is a statistical technique used to investigate the connection between a dependent variable and one or more independent variables.\n", - "\n", - "It is frequently employed in climatology to analyze trends and patterns in climatic data, identify correlations between different climatic parameters, and create models that can predict future changes. By identifying patterns and connections in massive datasets, regression analysis offers several benefits for weather research. For instance, regression analysis can be used to pinpoint the elements that affect global temperatures, such as solar radiation, atmospheric greenhouse gases, and volcanic eruptions. Climate scientists can create models that can accurately predict future changes by including these variables in a regression model.\n", - "\n", - "Moreover, regression analysis can assist climate experts in spotting natural fluctuations in climate data, like El Niño events, and in determining how human activities like deforestation and fossil fuel combustion affect the environment. Regression analysis can also evaluate the effectiveness of various mitigation tactics, such as carbon pricing policies or renewable energy initiatives.\n", - "\n", - "Overall, regression analysis is a potent tool for analyzing complex climate data and producing reliable projections of upcoming alterations.\n", - "\n", - ".. seealso::\n", - " - Regression Analysis: Definition, Types, Usage & Advantages. Website: https://www.questionpro.com/blog/regression-analysis/\n", - " - The Advantages of Regression Analysis & Forecasting. Website: https://smallbusiness.chron.com/advantages-regression-analysis-forecasting-61800.html\n", - "\n", - "In this subsection we try to regress the Niño 3.4 index on the Barents-Kara December SIC data.\n", - "Before performing the regression analysis, we can see that the longitude range of the SST data is **-180°~180°**,\n", - "try to convert the longitude range to **0°~360°** using :py:func:`easyclimate.utility.transfer_xarray_lon_from180TO360 `.\n", - "\n" + "## Regression\nRegression analysis is a statistical technique used to investigate the connection between a dependent variable and one or more independent variables.\n\nIt is frequently employed in climatology to analyze trends and patterns in climatic data, identify correlations between different climatic parameters, and create models that can predict future changes. By identifying patterns and connections in massive datasets, regression analysis offers several benefits for weather research. For instance, regression analysis can be used to pinpoint the elements that affect global temperatures, such as solar radiation, atmospheric greenhouse gases, and volcanic eruptions. Climate scientists can create models that can accurately predict future changes by including these variables in a regression model.\n\nMoreover, regression analysis can assist climate experts in spotting natural fluctuations in climate data, like El Ni\u00f1o events, and in determining how human activities like deforestation and fossil fuel combustion affect the environment. Regression analysis can also evaluate the effectiveness of various mitigation tactics, such as carbon pricing policies or renewable energy initiatives.\n\nOverall, regression analysis is a potent tool for analyzing complex climate data and producing reliable projections of upcoming alterations.\n\n.. seealso::\n - Regression Analysis: Definition, Types, Usage & Advantages. Website: https://www.questionpro.com/blog/regression-analysis/\n - The Advantages of Regression Analysis & Forecasting. Website: https://smallbusiness.chron.com/advantages-regression-analysis-forecasting-61800.html\n\nIn this subsection we try to regress the Ni\u00f1o 3.4 index on the Barents-Kara December SIC data.\nBefore performing the regression analysis, we can see that the longitude range of the SST data is **-180\u00b0~180\u00b0**,\ntry to convert the longitude range to **0\u00b0~360\u00b0** using :py:func:`easyclimate.utility.transfer_xarray_lon_from180TO360 `.\n\n" ] }, { @@ -296,21 +177,14 @@ }, "outputs": [], "source": [ - "sst_data_0_360 = ecl.utility.transfer_xarray_lon_from180TO360(sst_data)\n", - "sst_data_0_360" + "sst_data_0_360 = ecl.utility.transfer_xarray_lon_from180TO360(sst_data)\nsst_data_0_360" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Further, :py:func:`easyclimate.remove_seasonal_cycle_mean ` is used to remove the climate state of each\n", - "month in order to obtain the individual month anomalies.\n", - "The figure below illustrates the November SST anomaly in the tropical equatorial Pacific during the 1982-83 super El Niño.\n", - "\n", - ".. seealso::\n", - " Philander, S. Meteorology: Anomalous El Niño of 1982–83. Nature 305, 16 (1983). https://doi.org/10.1038/305016a0\n", - "\n" + "Further, :py:func:`easyclimate.remove_seasonal_cycle_mean ` is used to remove the climate state of each\nmonth in order to obtain the individual month anomalies.\nThe figure below illustrates the November SST anomaly in the tropical equatorial Pacific during the 1982-83 super El Ni\u00f1o.\n\n.. seealso::\n Philander, S. Meteorology: Anomalous El Ni\u00f1o of 1982\u201383. Nature 305, 16 (1983). https://doi.org/10.1038/305016a0\n\n" ] }, { @@ -321,34 +195,14 @@ }, "outputs": [], "source": [ - "sst_data_anormaly = ecl.remove_seasonal_cycle_mean(sst_data_0_360)\n", - "\n", - "fig, ax = plt.subplots(\n", - " figsize=(10, 4), subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "sst_data_anormaly.sel(lon=slice(120, 290)).isel(time=22).plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"bottom\", \"pad\": 0.1},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)" + "sst_data_anormaly = ecl.remove_seasonal_cycle_mean(sst_data_0_360)\n\nfig, ax = plt.subplots(\n figsize=(10, 4), subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nsst_data_anormaly.sel(lon=slice(120, 290)).isel(time=22).plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"bottom\", \"pad\": 0.1},\n cmap=\"RdBu_r\",\n levels=21,\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The Niño3.4 index is commonly used as an indicator for detecting ENSO,\n", - "and `easyclimate` provides :py:func:`easyclimate.field.air_sea_interaction.calc_index_nino34 ` to calculate the index using SST original dataset.\n", - "\n", - ".. seealso::\n", - " Anthony G. Bamston, Muthuvel Chelliah & Stanley B. Goldenberg (1997) Documentation of a highly ENSO‐related sst region in the equatorial pacific: Research note, Atmosphere-Ocean, 35:3, 367-383, DOI: https://doi.org/10.1080/07055900.1997.9649597\n", - "\n" + "The Ni\u00f1o3.4 index is commonly used as an indicator for detecting ENSO,\nand `easyclimate` provides :py:func:`easyclimate.field.air_sea_interaction.calc_index_nino34 ` to calculate the index using SST original dataset.\n\n.. seealso::\n Anthony G. Bamston, Muthuvel Chelliah & Stanley B. Goldenberg (1997) Documentation of a highly ENSO\u2010related sst region in the equatorial pacific: Research note, Atmosphere-Ocean, 35:3, 367-383, DOI: https://doi.org/10.1080/07055900.1997.9649597\n\n" ] }, { @@ -359,19 +213,14 @@ }, "outputs": [], "source": [ - "nino34_monthly_index = ecl.field.air_sea_interaction.calc_index_nino34(sst_data_0_360)\n", - "\n", - "nino34_monthly_index.plot(\n", - " figsize=(8, 3),\n", - ")" + "nino34_monthly_index = ecl.field.air_sea_interaction.calc_index_nino34(sst_data_0_360)\n\nnino34_monthly_index.plot(\n figsize=(8, 3),\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - ":py:func:`easyclimate.calc_yearly_climatological_mean ` is then used to solve for the annual average of the monthly index data\n", - "\n" + ":py:func:`easyclimate.calc_yearly_climatological_mean ` is then used to solve for the annual average of the monthly index data\n\n" ] }, { @@ -382,18 +231,14 @@ }, "outputs": [], "source": [ - "nino34_12_index = ecl.get_specific_months_data(nino34_monthly_index, 12)\n", - "nino34_dec_yearly_index = ecl.calc_yearly_climatological_mean(nino34_12_index)\n", - "nino34_dec_yearly_index" + "nino34_12_index = ecl.get_specific_months_data(nino34_monthly_index, 12)\nnino34_dec_yearly_index = ecl.calc_yearly_climatological_mean(nino34_12_index)\nnino34_dec_yearly_index" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Unlike solving for linear trend without passing in `x`, regression analysis must use the parameter `x` to pass in the object to be regressed.\n", - "Care must be taken to ensure that the `time` dimensions are identical.\n", - "\n" + "Unlike solving for linear trend without passing in `x`, regression analysis must use the parameter `x` to pass in the object to be regressed.\nCare must be taken to ensure that the `time` dimensions are identical.\n\n" ] }, { @@ -404,19 +249,14 @@ }, "outputs": [], "source": [ - "sic_reg_nino34 = ecl.calc_linregress_spatial(\n", - " sic_data_Barents_Sea_12, x=nino34_dec_yearly_index.data\n", - ")\n", - "sic_reg_nino34 = sic_reg_nino34.compute()\n", - "sic_reg_nino34" + "sic_reg_nino34 = ecl.calc_linregress_spatial(\n sic_data_Barents_Sea_12, x=nino34_dec_yearly_index.data\n)\nsic_reg_nino34 = sic_reg_nino34.compute()\nsic_reg_nino34" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Here is an attempt to plot the results of the regression analysis.\n", - "\n" + "Here is an attempt to plot the results of the regression analysis.\n\n" ] }, { @@ -427,40 +267,14 @@ }, "outputs": [], "source": [ - "draw_sic_slope = sic_reg_nino34.slope\n", - "draw_sic_pvalue = sic_reg_nino34.pvalue\n", - "\n", - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "draw_sic_slope.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "\n", - "ecl.plot.draw_significant_area_contourf(\n", - " draw_sic_pvalue, ax=ax, thresh=0.05, transform=ccrs.PlateCarree()\n", - ")" + "draw_sic_slope = sic_reg_nino34.slope\ndraw_sic_pvalue = sic_reg_nino34.pvalue\n\nfig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\ndraw_sic_slope.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"RdBu_r\",\n levels=21,\n)\n\necl.plot.draw_significant_area_contourf(\n draw_sic_pvalue, ax=ax, thresh=0.05, transform=ccrs.PlateCarree()\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Detrend\n", - "Sea ice area shows an approximately linear trend of decreasing due to global warming.\n", - "We remove the linear trend from SIC in order to study the variability of SIC itself.\n", - "In addition, here we explore the differences between the trend followed by regional averaging and regional averaging followed by detrending approaches.\n", - "\n" + "## Detrend\nSea ice area shows an approximately linear trend of decreasing due to global warming.\nWe remove the linear trend from SIC in order to study the variability of SIC itself.\nIn addition, here we explore the differences between the trend followed by regional averaging and regional averaging followed by detrending approaches.\n\n" ] }, { @@ -471,21 +285,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_spatial_mean = sic_data_Barents_Sea_12.mean(dim=(\"lat\", \"lon\"))\n", - "sic_data_Barents_Sea_12_spatial_detrendmean = ecl.calc_detrend_spatial(\n", - " sic_data_Barents_Sea_12, time_dim=\"time\"\n", - ").mean(dim=(\"lat\", \"lon\"))\n", - "sic_data_Barents_Sea_12_time_detrendmean = ecl.calc_detrend_spatial(\n", - " sic_data_Barents_Sea_12_spatial_mean, time_dim=\"time\"\n", - ")" + "sic_data_Barents_Sea_12_spatial_mean = sic_data_Barents_Sea_12.mean(dim=(\"lat\", \"lon\"))\nsic_data_Barents_Sea_12_spatial_detrendmean = ecl.calc_detrend_spatial(\n sic_data_Barents_Sea_12, time_dim=\"time\"\n).mean(dim=(\"lat\", \"lon\"))\nsic_data_Barents_Sea_12_time_detrendmean = ecl.calc_detrend_spatial(\n sic_data_Barents_Sea_12_spatial_mean, time_dim=\"time\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The results show that there is no significant difference between these two detrending methods to study the variability of SIC in the Barents-Kara Seas.\n", - "\n" + "The results show that there is no significant difference between these two detrending methods to study the variability of SIC in the Barents-Kara Seas.\n\n" ] }, { @@ -496,48 +303,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(2, 1, sharex=True)\n", - "\n", - "sic_data_Barents_Sea_12_spatial_mean.plot(ax=ax[0])\n", - "ax[0].set_xlabel(\"\")\n", - "ax[0].set_title(\"Original\")\n", - "\n", - "\n", - "sic_data_Barents_Sea_12_spatial_detrendmean.plot(\n", - " ax=ax[1], label=\"detrend -> spatial mean\"\n", - ")\n", - "sic_data_Barents_Sea_12_time_detrendmean.plot(\n", - " ax=ax[1], ls=\"--\", label=\"spatial mean -> detrend\"\n", - ")\n", - "ax[1].set_xlabel(\"\")\n", - "ax[1].set_title(\"Detrend\")\n", - "ax[1].legend()" + "fig, ax = plt.subplots(2, 1, sharex=True)\n\nsic_data_Barents_Sea_12_spatial_mean.plot(ax=ax[0])\nax[0].set_xlabel(\"\")\nax[0].set_title(\"Original\")\n\n\nsic_data_Barents_Sea_12_spatial_detrendmean.plot(\n ax=ax[1], label=\"detrend -> spatial mean\"\n)\nsic_data_Barents_Sea_12_time_detrendmean.plot(\n ax=ax[1], ls=\"--\", label=\"spatial mean -> detrend\"\n)\nax[1].set_xlabel(\"\")\nax[1].set_title(\"Detrend\")\nax[1].legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Weighted Spatial Data\n", - "When calculating regional averages in high-latitude areas,\n", - "considering weights is necessary because regions at different latitudes cover unequal\n", - "surface areas on the Earth. Since the Earth is approximately a spheroid, areas\n", - "closer to the poles have a different distribution of surface area on the spherical surface.\n", - "\n", - "One common way to incorporate weights is by using the cosine of latitude, i.e., multiplying by $\\cos (\\varphi)$,\n", - "where $\\varphi$ represents the latitude of a location. This is because areas at higher latitudes,\n", - "close to the poles, have higher latitudes and smaller cosine values, allowing for a\n", - "smaller weight to be applied to these regions when calculating averages.\n", - "\n", - "In summary, considering weights is done to more accurately account for the distribution\n", - "of surface area on the Earth, ensuring that contributions from different\n", - "regions are weighted according to their actual surface area when calculating\n", - "averages or other regional statistical measures.\n", - "\n", - ":py:func:`easyclimate.utility.get_weighted_spatial_data ` can help us create\n", - "an :py:class:`xarray.core.weighted.DataArrayWeighted ` object.\n", - "This object will automatically consider and calculate weights in subsequent area operations, thereby achieving the operation of the weighted spatial average.\n", - "\n" + "## Weighted Spatial Data\nWhen calculating regional averages in high-latitude areas,\nconsidering weights is necessary because regions at different latitudes cover unequal\nsurface areas on the Earth. Since the Earth is approximately a spheroid, areas\ncloser to the poles have a different distribution of surface area on the spherical surface.\n\nOne common way to incorporate weights is by using the cosine of latitude, i.e., multiplying by $\\cos (\\varphi)$,\nwhere $\\varphi$ represents the latitude of a location. This is because areas at higher latitudes,\nclose to the poles, have higher latitudes and smaller cosine values, allowing for a\nsmaller weight to be applied to these regions when calculating averages.\n\nIn summary, considering weights is done to more accurately account for the distribution\nof surface area on the Earth, ensuring that contributions from different\nregions are weighted according to their actual surface area when calculating\naverages or other regional statistical measures.\n\n:py:func:`easyclimate.utility.get_weighted_spatial_data ` can help us create\nan :py:class:`xarray.core.weighted.DataArrayWeighted ` object.\nThis object will automatically consider and calculate weights in subsequent area operations, thereby achieving the operation of the weighted spatial average.\n\n" ] }, { @@ -548,21 +321,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_detrend = ecl.calc_detrend_spatial(\n", - " sic_data_Barents_Sea_12, time_dim=\"time\"\n", - ")\n", - "grid_detrend_data_weighted_obj = ecl.utility.get_weighted_spatial_data(\n", - " sic_data_Barents_Sea_12_detrend, lat_dim=\"lat\", lon_dim=\"lon\"\n", - ")\n", - "print(type(grid_detrend_data_weighted_obj))" + "sic_data_Barents_Sea_12_detrend = ecl.calc_detrend_spatial(\n sic_data_Barents_Sea_12, time_dim=\"time\"\n)\ngrid_detrend_data_weighted_obj = ecl.utility.get_weighted_spatial_data(\n sic_data_Barents_Sea_12_detrend, lat_dim=\"lat\", lon_dim=\"lon\"\n)\nprint(type(grid_detrend_data_weighted_obj))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Solve for regional averaging for `grid_detrend_data_weighted_obj` objects (the role of weights is considered at this point)\n", - "\n" + "Solve for regional averaging for `grid_detrend_data_weighted_obj` objects (the role of weights is considered at this point)\n\n" ] }, { @@ -573,18 +339,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_spatial_detrend_weightedmean = (\n", - " grid_detrend_data_weighted_obj.mean(dim=(\"lat\", \"lon\"))\n", - ")\n", - "sic_data_Barents_Sea_12_spatial_detrend_weightedmean" + "sic_data_Barents_Sea_12_spatial_detrend_weightedmean = (\n grid_detrend_data_weighted_obj.mean(dim=(\"lat\", \"lon\"))\n)\nsic_data_Barents_Sea_12_spatial_detrend_weightedmean" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can find some differences between the data considering latitude weights and those not considering latitude weights.\n", - "\n" + "We can find some differences between the data considering latitude weights and those not considering latitude weights.\n\n" ] }, { @@ -595,50 +357,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(2, 1, sharex=True)\n", - "\n", - "sic_data_Barents_Sea_12_spatial_mean.plot(ax=ax[0])\n", - "ax[0].set_xlabel(\"\")\n", - "ax[0].set_title(\"Original\")\n", - "\n", - "\n", - "sic_data_Barents_Sea_12_spatial_detrendmean.plot(ax=ax[1], label=\"Regular mean\")\n", - "sic_data_Barents_Sea_12_spatial_detrend_weightedmean.plot(\n", - " ax=ax[1], ls=\"--\", label=\"Weighted mean\"\n", - ")\n", - "ax[1].set_xlabel(\"\")\n", - "ax[1].set_title(\"Detrend\")\n", - "ax[1].legend()\n", - "fig.show()" + "fig, ax = plt.subplots(2, 1, sharex=True)\n\nsic_data_Barents_Sea_12_spatial_mean.plot(ax=ax[0])\nax[0].set_xlabel(\"\")\nax[0].set_title(\"Original\")\n\n\nsic_data_Barents_Sea_12_spatial_detrendmean.plot(ax=ax[1], label=\"Regular mean\")\nsic_data_Barents_Sea_12_spatial_detrend_weightedmean.plot(\n ax=ax[1], ls=\"--\", label=\"Weighted mean\"\n)\nax[1].set_xlabel(\"\")\nax[1].set_title(\"Detrend\")\nax[1].legend()\nfig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Skewness\n", - "\n", - "Skewness is a measure of the asymmetry of a probability distribution.\n", - "\n", - "It quantifies the extent to which a distribution deviates from being symmetric around its mean. A distribution with a skewness\n", - "value of zero is considered symmetric, meaning that it has equal probabilities of occurring\n", - "above and below its mean. However, when the skewness value is non-zero, the distribution\n", - "becomes asymmetric, indicating that there is a higher likelihood of occurrence on one side\n", - "of the mean than the other.\n", - "\n", - "Skewness can be positive or negative, depending on whether the\n", - "distribution is skewed towards larger or smaller values.\n", - "\n", - "In climate analysis, skewness can arise due to various factors such as changes in atmospheric circulation patterns, uneven\n", - "temperature or precipitation distributions, or differences in measurement instruments.\n", - "\n", - ".. seealso::\n", - " - Distributions of Daily Meteorological Variables: Background. Website: https://psl.noaa.gov/data/atmoswrit/distributions/background/index.html\n", - " - Bakouch HS, Cadena M, Chesneau C. A new class of skew distributions with climate data analysis. J Appl Stat. 2020 Jul 13;48(16):3002-3024. doi: https://doi.org/10.1080/02664763.2020.1791804. PMID: 35707257; PMCID: PMC9042114.\n", - "\n", - "The skewness is calculated using :py:func:`easyclimate.calc_skewness_spatial `.\n", - "The result of the calculation contains the skewness and p-value.\n", - "\n" + "## Skewness\n\nSkewness is a measure of the asymmetry of a probability distribution.\n\nIt quantifies the extent to which a distribution deviates from being symmetric around its mean. A distribution with a skewness\nvalue of zero is considered symmetric, meaning that it has equal probabilities of occurring\nabove and below its mean. However, when the skewness value is non-zero, the distribution\nbecomes asymmetric, indicating that there is a higher likelihood of occurrence on one side\nof the mean than the other.\n\nSkewness can be positive or negative, depending on whether the\ndistribution is skewed towards larger or smaller values.\n\nIn climate analysis, skewness can arise due to various factors such as changes in atmospheric circulation patterns, uneven\ntemperature or precipitation distributions, or differences in measurement instruments.\n\n.. seealso::\n - Distributions of Daily Meteorological Variables: Background. Website: https://psl.noaa.gov/data/atmoswrit/distributions/background/index.html\n - Bakouch HS, Cadena M, Chesneau C. A new class of skew distributions with climate data analysis. J Appl Stat. 2020 Jul 13;48(16):3002-3024. doi: https://doi.org/10.1080/02664763.2020.1791804. PMID: 35707257; PMCID: PMC9042114.\n\nThe skewness is calculated using :py:func:`easyclimate.calc_skewness_spatial `.\nThe result of the calculation contains the skewness and p-value.\n\n" ] }, { @@ -649,13 +375,7 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_detrend = ecl.calc_detrend_spatial(\n", - " sic_data_Barents_Sea_12, time_dim=\"time\"\n", - ")\n", - "sic_data_Barents_Sea_12_skew = ecl.calc_skewness_spatial(\n", - " sic_data_Barents_Sea_12_detrend, dim=\"time\"\n", - ")\n", - "sic_data_Barents_Sea_12_skew" + "sic_data_Barents_Sea_12_detrend = ecl.calc_detrend_spatial(\n sic_data_Barents_Sea_12, time_dim=\"time\"\n)\nsic_data_Barents_Sea_12_skew = ecl.calc_skewness_spatial(\n sic_data_Barents_Sea_12_detrend, dim=\"time\"\n)\nsic_data_Barents_Sea_12_skew" ] }, { @@ -666,52 +386,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "# SIC slope\n", - "sic_data_Barents_Sea_12_skew.skewness.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "\n", - "# SIC 95% significant level\n", - "ecl.plot.draw_significant_area_contourf(\n", - " sic_data_Barents_Sea_12_skew.pvalue,\n", - " ax=ax,\n", - " thresh=0.05,\n", - " transform=ccrs.PlateCarree(),\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\n# SIC slope\nsic_data_Barents_Sea_12_skew.skewness.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"RdBu_r\",\n levels=21,\n)\n\n# SIC 95% significant level\necl.plot.draw_significant_area_contourf(\n sic_data_Barents_Sea_12_skew.pvalue,\n ax=ax,\n thresh=0.05,\n transform=ccrs.PlateCarree(),\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Kurtosis\n", - "Kurtosis is a measure of the \"tailedness\" of a probability distribution.\n", - "It describes how heavy or light the tails of a distribution are relative\n", - "to a standard normal distribution. A distribution with high kurtosis\n", - "has heavier tails and a greater propensity for extreme events, whereas a\n", - "distribution with low kurtosis has lighter tails and fewer extreme events.\n", - "Kurtosis is particularly useful in climate analysis because it can reveal\n", - "information about the frequency and intensity of extreme weather events such as hurricanes, droughts, or heatwaves.\n", - "\n", - ".. seealso::\n", - " - Distributions of Daily Meteorological Variables: Background. Website: https://psl.noaa.gov/data/atmoswrit/distributions/background/index.html\n", - " - Bakouch HS, Cadena M, Chesneau C. A new class of skew distributions with climate data analysis. J Appl Stat. 2020 Jul 13;48(16):3002-3024. doi: https://doi.org/10.1080/02664763.2020.1791804. PMID: 35707257; PMCID: PMC9042114.\n", - "\n", - "The skewness is calculated using :py:func:`easyclimate.calc_kurtosis_spatial `.\n", - "\n" + "## Kurtosis\nKurtosis is a measure of the \"tailedness\" of a probability distribution.\nIt describes how heavy or light the tails of a distribution are relative\nto a standard normal distribution. A distribution with high kurtosis\nhas heavier tails and a greater propensity for extreme events, whereas a\ndistribution with low kurtosis has lighter tails and fewer extreme events.\nKurtosis is particularly useful in climate analysis because it can reveal\ninformation about the frequency and intensity of extreme weather events such as hurricanes, droughts, or heatwaves.\n\n.. seealso::\n - Distributions of Daily Meteorological Variables: Background. Website: https://psl.noaa.gov/data/atmoswrit/distributions/background/index.html\n - Bakouch HS, Cadena M, Chesneau C. A new class of skew distributions with climate data analysis. J Appl Stat. 2020 Jul 13;48(16):3002-3024. doi: https://doi.org/10.1080/02664763.2020.1791804. PMID: 35707257; PMCID: PMC9042114.\n\nThe skewness is calculated using :py:func:`easyclimate.calc_kurtosis_spatial `.\n\n" ] }, { @@ -722,18 +404,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_kurt = ecl.calc_kurtosis_spatial(\n", - " sic_data_Barents_Sea_12_detrend, dim=\"time\"\n", - ")\n", - "sic_data_Barents_Sea_12_kurt" + "sic_data_Barents_Sea_12_kurt = ecl.calc_kurtosis_spatial(\n sic_data_Barents_Sea_12_detrend, dim=\"time\"\n)\nsic_data_Barents_Sea_12_kurt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Consider plotting kurtosis below\n", - "\n" + "Consider plotting kurtosis below\n\n" ] }, { @@ -744,51 +422,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "# SIC slope\n", - "sic_data_Barents_Sea_12_kurt.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\n# SIC slope\nsic_data_Barents_Sea_12_kurt.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"RdBu_r\",\n levels=21,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Composite Analysis\n", - "In the process of climate analysis, **composite analysis** is a statistical and integrative method\n", - "used to study the spatial and temporal distribution of specific climate events or phenomena.\n", - "The primary purpose of this analysis method is to identify common features and patterns\n", - "among multiple events or time periods by combining their information.\n", - "\n", - "Specifically, the steps of composite analysis typically include the following aspects:\n", - "\n", - "1. **Event Selection**: Firstly, a set of events related to the research objective is chosen. These events could be specific climate phenomena such as heavy rainfall, drought, or temperature anomalies.\n", - "2. **Data Collection**: Collect meteorological data related to the selected events, including observational data, model outputs, or remote sensing data.\n", - "3. **Event Alignment**: Time-align the chosen events to ensure that they are within the same temporal framework for analysis.\n", - "4. **Data Combination**: Combine data values at corresponding time points into a composite dataset. Averages or weighted averages are often used to reduce the impact of random noise.\n", - "\n", - "The advantages of this method include:\n", - "\n", - "- **Highlighting Common Features**: By combining data from multiple events or time periods, composite analysis can highlight common features, aiding in the identification of general patterns in climate events.\n", - "- **Noise Reduction**: By averaging data, composite analysis helps to reduce the impact of random noise, resulting in more stable and reliable analysis outcomes.\n", - "- **Spatial Consistency**: Through spatial averaging, this method helps reveal the consistent spatial distribution of climate events, providing a more comprehensive understanding.\n", - "- **Facilitating Comparisons**: Composite analysis makes it convenient to compare different events or time periods as it integrates them into a unified framework.\n", - "\n", - "Here we try to extract the El Niño and La Niña events using the standard deviation of the Niño 3.4 index as a threshold.\n", - "\n" + "## Composite Analysis\nIn the process of climate analysis, **composite analysis** is a statistical and integrative method\nused to study the spatial and temporal distribution of specific climate events or phenomena.\nThe primary purpose of this analysis method is to identify common features and patterns\namong multiple events or time periods by combining their information.\n\nSpecifically, the steps of composite analysis typically include the following aspects:\n\n1. **Event Selection**: Firstly, a set of events related to the research objective is chosen. These events could be specific climate phenomena such as heavy rainfall, drought, or temperature anomalies.\n2. **Data Collection**: Collect meteorological data related to the selected events, including observational data, model outputs, or remote sensing data.\n3. **Event Alignment**: Time-align the chosen events to ensure that they are within the same temporal framework for analysis.\n4. **Data Combination**: Combine data values at corresponding time points into a composite dataset. Averages or weighted averages are often used to reduce the impact of random noise.\n\nThe advantages of this method include:\n\n- **Highlighting Common Features**: By combining data from multiple events or time periods, composite analysis can highlight common features, aiding in the identification of general patterns in climate events.\n- **Noise Reduction**: By averaging data, composite analysis helps to reduce the impact of random noise, resulting in more stable and reliable analysis outcomes.\n- **Spatial Consistency**: Through spatial averaging, this method helps reveal the consistent spatial distribution of climate events, providing a more comprehensive understanding.\n- **Facilitating Comparisons**: Composite analysis makes it convenient to compare different events or time periods as it integrates them into a unified framework.\n\nHere we try to extract the El Ni\u00f1o and La Ni\u00f1a events using the standard deviation of the Ni\u00f1o 3.4 index as a threshold.\n\n" ] }, { @@ -799,17 +440,14 @@ }, "outputs": [], "source": [ - "nino34_dec_yearly_index_std = nino34_dec_yearly_index.std(dim=\"time\").data\n", - "nino34_dec_yearly_index_std" + "nino34_dec_yearly_index_std = nino34_dec_yearly_index.std(dim=\"time\").data\nnino34_dec_yearly_index_std" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - ":py:func:`easyclimate.get_year_exceed_index_upper_bound ` is able to obtain the years that exceed the upper bound of the Niño 3.4 exponential threshold,\n", - "and similarly :py:func:`easyclimate.get_year_exceed_index_lower_bound ` can obtain the years that exceed the lower bound of the exponential threshold.\n", - "\n" + ":py:func:`easyclimate.get_year_exceed_index_upper_bound ` is able to obtain the years that exceed the upper bound of the Ni\u00f1o 3.4 exponential threshold,\nand similarly :py:func:`easyclimate.get_year_exceed_index_lower_bound ` can obtain the years that exceed the lower bound of the exponential threshold.\n\n" ] }, { @@ -820,23 +458,14 @@ }, "outputs": [], "source": [ - "elnino_year = ecl.get_year_exceed_index_upper_bound(\n", - " nino34_dec_yearly_index, thresh=nino34_dec_yearly_index_std\n", - ")\n", - "lanina_year = ecl.get_year_exceed_index_lower_bound(\n", - " nino34_dec_yearly_index, thresh=-nino34_dec_yearly_index_std\n", - ")\n", - "print(\"El niño years: \", elnino_year)\n", - "print(\"La niña years: \", lanina_year)" + "elnino_year = ecl.get_year_exceed_index_upper_bound(\n nino34_dec_yearly_index, thresh=nino34_dec_yearly_index_std\n)\nlanina_year = ecl.get_year_exceed_index_lower_bound(\n nino34_dec_yearly_index, thresh=-nino34_dec_yearly_index_std\n)\nprint(\"El ni\u00f1o years: \", elnino_year)\nprint(\"La ni\u00f1a years: \", lanina_year)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Further we use :py:func:`easyclimate.get_specific_years_data ` to extract data\n", - "for the El Niño years within `sic_data_Barents_Sea_12_detrend`. The results show that six temporal levels were extracted.\n", - "\n" + "Further we use :py:func:`easyclimate.get_specific_years_data ` to extract data\nfor the El Ni\u00f1o years within `sic_data_Barents_Sea_12_detrend`. The results show that six temporal levels were extracted.\n\n" ] }, { @@ -847,23 +476,14 @@ }, "outputs": [], "source": [ - "sic_data_Barents_Sea_12_detrend_elnino = ecl.get_specific_years_data(\n", - " sic_data_Barents_Sea_12_detrend, elnino_year\n", - ")\n", - "sic_data_Barents_Sea_12_detrend_elnino.name = \"El niño years\"\n", - "sic_data_Barents_Sea_12_detrend_lanina = ecl.get_specific_years_data(\n", - " sic_data_Barents_Sea_12_detrend, lanina_year\n", - ")\n", - "sic_data_Barents_Sea_12_detrend_lanina.name = \"La niña years\"\n", - "sic_data_Barents_Sea_12_detrend_elnino" + "sic_data_Barents_Sea_12_detrend_elnino = ecl.get_specific_years_data(\n sic_data_Barents_Sea_12_detrend, elnino_year\n)\nsic_data_Barents_Sea_12_detrend_elnino.name = \"El ni\u00f1o years\"\nsic_data_Barents_Sea_12_detrend_lanina = ecl.get_specific_years_data(\n sic_data_Barents_Sea_12_detrend, lanina_year\n)\nsic_data_Barents_Sea_12_detrend_lanina.name = \"La ni\u00f1a years\"\nsic_data_Barents_Sea_12_detrend_elnino" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Similarly, 5 temporal levels were extracted for the La Niña years.\n", - "\n" + "Similarly, 5 temporal levels were extracted for the La Ni\u00f1a years.\n\n" ] }, { @@ -881,8 +501,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We now plot the distribution of SIC during El Niño and La Niña years.\n", - "\n" + "We now plot the distribution of SIC during El Ni\u00f1o and La Ni\u00f1a years.\n\n" ] }, { @@ -893,48 +512,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " 1,\n", - " 2,\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " },\n", - " figsize=(10, 4),\n", - ")\n", - "\n", - "for axi in ax.flat:\n", - " axi.gridlines(\n", - " draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\"\n", - " )\n", - " axi.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "sic_data_Barents_Sea_12_detrend_elnino.mean(dim=\"time\").plot.contourf(\n", - " ax=ax[0],\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"bottom\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "ax[0].set_title(\"El niño years\")\n", - "\n", - "sic_data_Barents_Sea_12_detrend_lanina.mean(dim=\"time\").plot.contourf(\n", - " ax=ax[1],\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"bottom\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "ax[1].set_title(\"La niña years\")" + "fig, ax = plt.subplots(\n 1,\n 2,\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n },\n figsize=(10, 4),\n)\n\nfor axi in ax.flat:\n axi.gridlines(\n draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\"\n )\n axi.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nsic_data_Barents_Sea_12_detrend_elnino.mean(dim=\"time\").plot.contourf(\n ax=ax[0],\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"bottom\"},\n cmap=\"RdBu_r\",\n levels=21,\n)\nax[0].set_title(\"El ni\u00f1o years\")\n\nsic_data_Barents_Sea_12_detrend_lanina.mean(dim=\"time\").plot.contourf(\n ax=ax[1],\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"bottom\"},\n cmap=\"RdBu_r\",\n levels=21,\n)\nax[1].set_title(\"La ni\u00f1a years\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "So is there a significant difference in the distribution of SIC between these two events?\n", - ":py:func:`easyclimate.calc_ttestSpatialPattern_spatial ` provides a\n", - "two-sample t-test operation to investigate whether there is a significant difference between the means of the two samples.\n", - "\n" + "So is there a significant difference in the distribution of SIC between these two events?\n:py:func:`easyclimate.calc_ttestSpatialPattern_spatial ` provides a\ntwo-sample t-test operation to investigate whether there is a significant difference between the means of the two samples.\n\n" ] }, { @@ -945,20 +530,14 @@ }, "outputs": [], "source": [ - "sig_diff = ecl.calc_ttestSpatialPattern_spatial(\n", - " sic_data_Barents_Sea_12_detrend_elnino,\n", - " sic_data_Barents_Sea_12_detrend_lanina,\n", - " dim=\"time\",\n", - ")\n", - "sig_diff" + "sig_diff = ecl.calc_ttestSpatialPattern_spatial(\n sic_data_Barents_Sea_12_detrend_elnino,\n sic_data_Barents_Sea_12_detrend_lanina,\n dim=\"time\",\n)\nsig_diff" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can find that there is little difference in the effect on SIC under different ENSO events.\n", - "\n" + "We can find that there is little difference in the effect on SIC under different ENSO events.\n\n" ] }, { @@ -969,32 +548,7 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\n", - " \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n", - " }\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "# SIC slope\n", - "diff = sic_data_Barents_Sea_12_detrend_lanina.mean(\n", - " dim=\"time\"\n", - ") - sic_data_Barents_Sea_12_detrend_elnino.mean(dim=\"time\")\n", - "diff.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " cbar_kwargs={\"location\": \"right\"},\n", - " cmap=\"RdBu_r\",\n", - " levels=21,\n", - ")\n", - "\n", - "ax.set_title(\"La niña minus El niño\", loc=\"left\")\n", - "\n", - "ecl.plot.draw_significant_area_contourf(\n", - " sig_diff.pvalue, ax=ax, thresh=0.1, transform=ccrs.PlateCarree()\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\n \"projection\": ccrs.Orthographic(central_longitude=70, central_latitude=70)\n }\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\n# SIC slope\ndiff = sic_data_Barents_Sea_12_detrend_lanina.mean(\n dim=\"time\"\n) - sic_data_Barents_Sea_12_detrend_elnino.mean(dim=\"time\")\ndiff.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n cbar_kwargs={\"location\": \"right\"},\n cmap=\"RdBu_r\",\n levels=21,\n)\n\nax.set_title(\"La ni\u00f1a minus El ni\u00f1o\", loc=\"left\")\n\necl.plot.draw_significant_area_contourf(\n sig_diff.pvalue, ax=ax, thresh=0.1, transform=ccrs.PlateCarree()\n)" ] } ], @@ -1014,7 +568,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/docs/example/plot_formatting_coordinates.ipynb b/docs/example/plot_formatting_coordinates.ipynb index 83d57d7c..dcb88c6c 100644 --- a/docs/example/plot_formatting_coordinates.ipynb +++ b/docs/example/plot_formatting_coordinates.ipynb @@ -251,7 +251,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/docs/example/plot_geographic_finite_difference.ipynb b/docs/example/plot_geographic_finite_difference.ipynb index a5728a12..bfac3a87 100644 --- a/docs/example/plot_geographic_finite_difference.ipynb +++ b/docs/example/plot_geographic_finite_difference.ipynb @@ -4,10 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "\n", - "# Geographic Finite Difference\n", - "\n", - "Before proceeding with all the steps, first import some necessary libraries and packages\n" + "\n# Geographic Finite Difference\n\nBefore proceeding with all the steps, first import some necessary libraries and packages\n" ] }, { @@ -18,18 +15,14 @@ }, "outputs": [], "source": [ - "import easyclimate as ecl\n", - "import xarray as xr\n", - "import matplotlib.pyplot as plt\n", - "import cartopy.crs as ccrs" + "import easyclimate as ecl\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Then consider obtaining meridional and zonal wind variables in tutorial data\n", - "\n" + "Then consider obtaining meridional and zonal wind variables in tutorial data\n\n" ] }, { @@ -40,29 +33,14 @@ }, "outputs": [], "source": [ - "u_data = ecl.tutorial.open_tutorial_dataset(\"uwnd_202201_mon_mean\").sortby(\"lat\").uwnd\n", - "v_data = ecl.tutorial.open_tutorial_dataset(\"vwnd_202201_mon_mean\").sortby(\"lat\").vwnd\n", - "z_data = ecl.tutorial.open_tutorial_dataset(\"hgt_202201_mon_mean\").sortby(\"lat\").hgt\n", - "temp_data = ecl.tutorial.open_tutorial_dataset(\"air_202201_mon_mean\").sortby(\"lat\").air\n", - "q_data = ecl.tutorial.open_tutorial_dataset(\"shum_202201_mon_mean\").sortby(\"lat\").shum\n", - "msl_data = (\n", - " ecl.tutorial.open_tutorial_dataset(\"pressfc_202201_mon_mean\").sortby(\"lat\").pres\n", - ")\n", - "pr_data = (\n", - " ecl.tutorial.open_tutorial_dataset(\"precip_202201_mon_mean\").sortby(\"lat\").precip\n", - ")\n", - "\n", - "uvdata = xr.Dataset()\n", - "uvdata[\"uwnd\"] = u_data\n", - "uvdata[\"vwnd\"] = v_data" + "u_data = ecl.tutorial.open_tutorial_dataset(\"uwnd_202201_mon_mean\").sortby(\"lat\").uwnd\nv_data = ecl.tutorial.open_tutorial_dataset(\"vwnd_202201_mon_mean\").sortby(\"lat\").vwnd\nz_data = ecl.tutorial.open_tutorial_dataset(\"hgt_202201_mon_mean\").sortby(\"lat\").hgt\ntemp_data = ecl.tutorial.open_tutorial_dataset(\"air_202201_mon_mean\").sortby(\"lat\").air\nq_data = ecl.tutorial.open_tutorial_dataset(\"shum_202201_mon_mean\").sortby(\"lat\").shum\nmsl_data = (\n ecl.tutorial.open_tutorial_dataset(\"pressfc_202201_mon_mean\").sortby(\"lat\").pres\n)\npr_data = (\n ecl.tutorial.open_tutorial_dataset(\"precip_202201_mon_mean\").sortby(\"lat\").precip\n)\n\nuvdata = xr.Dataset()\nuvdata[\"uwnd\"] = u_data\nuvdata[\"vwnd\"] = v_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Obtain data slices on 500hPa isobars for January 2022\n", - "\n" + "Obtain data slices on 500hPa isobars for January 2022\n\n" ] }, { @@ -73,17 +51,14 @@ }, "outputs": [], "source": [ - "uvdata_500_202201 = uvdata.sel(level=500, time=\"2022-01-01\")\n", - "z_data_500_202201 = z_data.sel(level=500, time=\"2022-01-01\")\n", - "temp_data_500_202201 = temp_data.sel(level=500, time=\"2022-01-01\")" + "uvdata_500_202201 = uvdata.sel(level=500, time=\"2022-01-01\")\nz_data_500_202201 = z_data.sel(level=500, time=\"2022-01-01\")\ntemp_data_500_202201 = temp_data.sel(level=500, time=\"2022-01-01\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Plotting a sample `quiver` plot of this data slice\n", - "\n" + "Plotting a sample `quiver` plot of this data slice\n\n" ] }, { @@ -94,37 +69,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.stock_img()\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "uvdata_500_202201.thin(lon=3, lat=3).plot.quiver(\n", - " ax=ax,\n", - " u=\"uwnd\",\n", - " v=\"vwnd\",\n", - " x=\"lon\",\n", - " y=\"lat\",\n", - " # projection on data\n", - " transform=ccrs.PlateCarree(),\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.stock_img()\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nuvdata_500_202201.thin(lon=3, lat=3).plot.quiver(\n ax=ax,\n u=\"uwnd\",\n v=\"vwnd\",\n x=\"lon\",\n y=\"lat\",\n # projection on data\n transform=ccrs.PlateCarree(),\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## First-order Partial Derivative\n", - "\n", - "Consider the function :py:func:`easyclimate.calc_gradient ` to compute the gradient of the zonal wind with respect to longitude.\n", - "\n", - "\\begin{align}\\frac{\\partial u}{\\partial \\lambda}\\end{align}\n", - "\n", - "The argument `dim` to the function :py:func:`easyclimate.calc_gradient ` specifies that the direction of the solution is `longitude`.\n", - "\n" + "## First-order Partial Derivative\n\nConsider the function :py:func:`easyclimate.calc_gradient ` to compute the gradient of the zonal wind with respect to longitude.\n\n\\begin{align}\\frac{\\partial u}{\\partial \\lambda}\\end{align}\n\nThe argument `dim` to the function :py:func:`easyclimate.calc_gradient ` specifies that the direction of the solution is `longitude`.\n\n" ] }, { @@ -135,9 +87,7 @@ }, "outputs": [], "source": [ - "uwnd_dx = ecl.calc_gradient(uvdata_500_202201.uwnd, dim=\"lon\")\n", - "\n", - "uwnd_dx" + "uwnd_dx = ecl.calc_gradient(uvdata_500_202201.uwnd, dim=\"lon\")\n\nuwnd_dx" ] }, { @@ -148,30 +98,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.stock_img()\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "uwnd_dx.plot.contourf(\n", - " ax=ax,\n", - " # projection on data\n", - " transform=ccrs.PlateCarree(),\n", - " # Colorbar is placed at the bottom\n", - " cbar_kwargs={\"location\": \"bottom\"},\n", - " levels=21,\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.stock_img()\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nuwnd_dx.plot.contourf(\n ax=ax,\n # projection on data\n transform=ccrs.PlateCarree(),\n # Colorbar is placed at the bottom\n cbar_kwargs={\"location\": \"bottom\"},\n levels=21,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Of course, it is also possible to pass in :py:class:`xarray.Dataset` directly into the function :py:func:`easyclimate.calc_gradient ` to iterate through all the variables, so that you can get the gradient of both the zonal and meridional winds with respect to longitude at the same time.\n", - "\n" + "Of course, it is also possible to pass in :py:class:`xarray.Dataset` directly into the function :py:func:`easyclimate.calc_gradient ` to iterate through all the variables, so that you can get the gradient of both the zonal and meridional winds with respect to longitude at the same time.\n\n" ] }, { @@ -182,19 +116,14 @@ }, "outputs": [], "source": [ - "uvwnd_dx = ecl.calc_gradient(uvdata_500_202201, dim=\"lon\")\n", - "\n", - "uvwnd_dx" + "uvwnd_dx = ecl.calc_gradient(uvdata_500_202201, dim=\"lon\")\n\nuvwnd_dx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "However, if one is required to solve for the gradient of the zonal wind with respect to the corresponding distance at each longitude, the function `calc_lon_gradient` should be used to calculate.\n", - "\n", - "\\begin{align}\\frac{\\partial F}{\\partial x} = \\frac{1}{R \\cos\\varphi} \\cdot \\frac{\\partial F}{\\partial \\lambda}\\end{align}\n", - "\n" + "However, if one is required to solve for the gradient of the zonal wind with respect to the corresponding distance at each longitude, the function `calc_lon_gradient` should be used to calculate.\n\n\\begin{align}\\frac{\\partial F}{\\partial x} = \\frac{1}{R \\cos\\varphi} \\cdot \\frac{\\partial F}{\\partial \\lambda}\\end{align}\n\n" ] }, { @@ -205,45 +134,21 @@ }, "outputs": [], "source": [ - "uwnd_dlon = ecl.calc_lon_gradient(uvdata_500_202201.uwnd, lon_dim=\"lon\", lat_dim=\"lat\")\n", - "\n", - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "uwnd_dlon.plot.contourf(\n", - " ax=ax,\n", - " # projection on data\n", - " transform=ccrs.PlateCarree(),\n", - " # Colorbar is placed at the bottom\n", - " cbar_kwargs={\"location\": \"bottom\"},\n", - " levels=21,\n", - ")" + "uwnd_dlon = ecl.calc_lon_gradient(uvdata_500_202201.uwnd, lon_dim=\"lon\", lat_dim=\"lat\")\n\nfig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nuwnd_dlon.plot.contourf(\n ax=ax,\n # projection on data\n transform=ccrs.PlateCarree(),\n # Colorbar is placed at the bottom\n cbar_kwargs={\"location\": \"bottom\"},\n levels=21,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Similarly, use :py:func:`easyclimate.calc_lat_gradient ` to solve for the gradient of the meridional wind with respect to the corresponding distance at each latitude.\n", - "\n" + "Similarly, use :py:func:`easyclimate.calc_lat_gradient ` to solve for the gradient of the meridional wind with respect to the corresponding distance at each latitude.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Second-order Partial Derivative\n", - "\n", - "The solution of the second-order partial derivative relies on three functional calculations\n", - "\n", - "- :py:func:`easyclimate.calc_lon_laplacian `: calculation of the second-order partial derivative term (Laplace term) along longitude.\n", - "\n", - "\\begin{align}\\frac{\\partial^2 F}{\\partial x^2} = \\frac{1}{(R \\cos\\varphi)^2} \\cdot \\frac{\\partial^2 F}{\\partial \\lambda^2}\\end{align}\n", - "\n" + "## Second-order Partial Derivative\n\nThe solution of the second-order partial derivative relies on three functional calculations\n\n- :py:func:`easyclimate.calc_lon_laplacian `: calculation of the second-order partial derivative term (Laplace term) along longitude.\n\n\\begin{align}\\frac{\\partial^2 F}{\\partial x^2} = \\frac{1}{(R \\cos\\varphi)^2} \\cdot \\frac{\\partial^2 F}{\\partial \\lambda^2}\\end{align}\n\n" ] }, { @@ -254,19 +159,14 @@ }, "outputs": [], "source": [ - "uwnd_dlon2 = ecl.calc_lon_laplacian(\n", - " uvdata_500_202201.uwnd, lon_dim=\"lon\", lat_dim=\"lat\"\n", - ")" + "uwnd_dlon2 = ecl.calc_lon_laplacian(\n uvdata_500_202201.uwnd, lon_dim=\"lon\", lat_dim=\"lat\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "- :py:func:`easyclimate.calc_lat_laplacian `: calculation of the second-order partial derivative term (Laplace term) along latitude.\n", - "\n", - "\\begin{align}\\frac{\\partial^2 F}{\\partial y^2} = \\frac{1}{R^2} \\cdot \\frac{\\partial^2 F}{\\partial \\varphi^2}\\end{align}\n", - "\n" + "- :py:func:`easyclimate.calc_lat_laplacian `: calculation of the second-order partial derivative term (Laplace term) along latitude.\n\n\\begin{align}\\frac{\\partial^2 F}{\\partial y^2} = \\frac{1}{R^2} \\cdot \\frac{\\partial^2 F}{\\partial \\varphi^2}\\end{align}\n\n" ] }, { @@ -284,11 +184,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "- :py:func:`easyclimate.calc_lon_lat_mixed_derivatives `: second-order mixed partial derivative terms along longitude and latitude.\n", - "\n", - "\\begin{align}\\frac{\\partial^2 F}{\\partial x \\partial y} = \\frac{1}{R^2 \\cos\\varphi} \\cdot \\frac{\\partial^2 F}{\\partial \\lambda \\partial \\varphi}\\end{align}\n", - "\n", - "\n" + "- :py:func:`easyclimate.calc_lon_lat_mixed_derivatives `: second-order mixed partial derivative terms along longitude and latitude.\n\n\\begin{align}\\frac{\\partial^2 F}{\\partial x \\partial y} = \\frac{1}{R^2 \\cos\\varphi} \\cdot \\frac{\\partial^2 F}{\\partial \\lambda \\partial \\varphi}\\end{align}\n\n\n" ] }, { @@ -299,17 +195,14 @@ }, "outputs": [], "source": [ - "uwnd_dlonlat = ecl.calc_lon_lat_mixed_derivatives(\n", - " uvdata_500_202201.uwnd, lon_dim=\"lon\", lat_dim=\"lat\"\n", - ")" + "uwnd_dlonlat = ecl.calc_lon_lat_mixed_derivatives(\n uvdata_500_202201.uwnd, lon_dim=\"lon\", lat_dim=\"lat\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Second-order partial derivative term along longitude.\n", - "\n" + "Second-order partial derivative term along longitude.\n\n" ] }, { @@ -320,25 +213,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "uwnd_dlon2.plot.contourf(\n", - " ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n", - ")\n", - "ax.set_title(\"$\\\\frac{\\\\partial^2 F}{\\\\partial x^2}$\", fontsize=20)" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nuwnd_dlon2.plot.contourf(\n ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n)\nax.set_title(\"$\\\\frac{\\\\partial^2 F}{\\\\partial x^2}$\", fontsize=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Second-order partial derivative term along latitude.\n", - "\n" + "Second-order partial derivative term along latitude.\n\n" ] }, { @@ -349,25 +231,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "uwnd_dlat2.plot.contourf(\n", - " ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n", - ")\n", - "ax.set_title(\"$\\\\frac{\\\\partial^2 F}{\\\\partial y^2}$\", fontsize=20)" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nuwnd_dlat2.plot.contourf(\n ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n)\nax.set_title(\"$\\\\frac{\\\\partial^2 F}{\\\\partial y^2}$\", fontsize=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Second-order mixed partial derivative terms along longitude and latitude.\n", - "\n" + "Second-order mixed partial derivative terms along longitude and latitude.\n\n" ] }, { @@ -378,32 +249,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "uwnd_dlonlat.plot.contourf(\n", - " ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n", - ")\n", - "ax.set_title(\"$\\\\frac{\\\\partial^2 F}{\\\\partial x \\\\partial y}$\", fontsize=20)" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\nuwnd_dlonlat.plot.contourf(\n ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n)\nax.set_title(\"$\\\\frac{\\\\partial^2 F}{\\\\partial x \\\\partial y}$\", fontsize=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Vorticity and Divergence\n", - "\n", - "Vorticity and divergence are measures of the degree of atmospheric rotation and volumetric flux per unit volume respectively. For vorticity and divergence in the quasi-geostrophic case, the potential height is used as input data for the calculations. In general, we first calculate the quasi-geostrophic wind.\n", - "\n", - "- :py:func:`easyclimate.calc_geostrophic_wind `: calculate the geostrophic wind.\n", - "\n", - "\\begin{align}u_g = - \\frac{g}{f} \\frac{\\partial H}{\\partial y}, \\ v_g = \\frac{g}{f} \\frac{\\partial H}{\\partial x}\\end{align}\n", - "\n", - "\n" + "## Vorticity and Divergence\n\nVorticity and divergence are measures of the degree of atmospheric rotation and volumetric flux per unit volume respectively. For vorticity and divergence in the quasi-geostrophic case, the potential height is used as input data for the calculations. In general, we first calculate the quasi-geostrophic wind.\n\n- :py:func:`easyclimate.calc_geostrophic_wind `: calculate the geostrophic wind.\n\n\\begin{align}u_g = - \\frac{g}{f} \\frac{\\partial H}{\\partial y}, \\ v_g = \\frac{g}{f} \\frac{\\partial H}{\\partial x}\\end{align}\n\n\n" ] }, { @@ -414,22 +267,14 @@ }, "outputs": [], "source": [ - "geostrophic_wind_data_500_202201 = ecl.calc_geostrophic_wind(\n", - " z_data_500_202201, lon_dim=\"lon\", lat_dim=\"lat\"\n", - ")" + "geostrophic_wind_data_500_202201 = ecl.calc_geostrophic_wind(\n z_data_500_202201, lon_dim=\"lon\", lat_dim=\"lat\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The function :py:func:`easyclimate.calc_vorticity ` is then used to compute the quasi-geostrophic vorticity.\n", - "\n", - "- :py:func:`easyclimate.calc_vorticity `: calculate the horizontal relative vorticity term.\n", - "\n", - "\\begin{align}\\zeta = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u}{R} \\tan \\varphi\\end{align}\n", - "\n", - "\n" + "The function :py:func:`easyclimate.calc_vorticity ` is then used to compute the quasi-geostrophic vorticity.\n\n- :py:func:`easyclimate.calc_vorticity `: calculate the horizontal relative vorticity term.\n\n\\begin{align}\\zeta = \\frac{\\partial v}{\\partial x} - \\frac{\\partial u}{\\partial y} + \\frac{u}{R} \\tan \\varphi\\end{align}\n\n\n" ] }, { @@ -440,22 +285,14 @@ }, "outputs": [], "source": [ - "qg_vor_data_500_202201 = ecl.calc_vorticity(\n", - " u_data=geostrophic_wind_data_500_202201.ug,\n", - " v_data=geostrophic_wind_data_500_202201.vg,\n", - " lon_dim=\"lon\",\n", - " lat_dim=\"lat\",\n", - ")\n", - "\n", - "qg_vor_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" + "qg_vor_data_500_202201 = ecl.calc_vorticity(\n u_data=geostrophic_wind_data_500_202201.ug,\n v_data=geostrophic_wind_data_500_202201.vg,\n lon_dim=\"lon\",\n lat_dim=\"lat\",\n)\n\nqg_vor_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Similar vorticity for actual winds, but for actual winds rather than quasi-geostrophic winds.\n", - "\n" + "Similar vorticity for actual winds, but for actual winds rather than quasi-geostrophic winds.\n\n" ] }, { @@ -466,28 +303,14 @@ }, "outputs": [], "source": [ - "vor_data_500_202201 = ecl.calc_vorticity(\n", - " u_data=uvdata_500_202201[\"uwnd\"],\n", - " v_data=uvdata_500_202201[\"vwnd\"],\n", - " lon_dim=\"lon\",\n", - " lat_dim=\"lat\",\n", - ")\n", - "\n", - "vor_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" + "vor_data_500_202201 = ecl.calc_vorticity(\n u_data=uvdata_500_202201[\"uwnd\"],\n v_data=uvdata_500_202201[\"vwnd\"],\n lon_dim=\"lon\",\n lat_dim=\"lat\",\n)\n\nvor_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In addition, the function :py:func:`easyclimate.calc_divergence ` calculate the quasi-geostrophic divergence.\n", - "\n", - "\\begin{align}\\mathrm{D} = \\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y} - \\frac{v}{R} \\tan \\varphi\\end{align}\n", - "\n", - "- :py:func:`easyclimate.calc_divergence `: calculate the horizontal divergence term.\n", - "\n", - "Quasi-geostrophic divergence\n", - "\n" + "In addition, the function :py:func:`easyclimate.calc_divergence ` calculate the quasi-geostrophic divergence.\n\n\\begin{align}\\mathrm{D} = \\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y} - \\frac{v}{R} \\tan \\varphi\\end{align}\n\n- :py:func:`easyclimate.calc_divergence `: calculate the horizontal divergence term.\n\nQuasi-geostrophic divergence\n\n" ] }, { @@ -498,22 +321,14 @@ }, "outputs": [], "source": [ - "qg_div_data_500_202201 = ecl.calc_divergence(\n", - " u_data=geostrophic_wind_data_500_202201.ug,\n", - " v_data=geostrophic_wind_data_500_202201.vg,\n", - " lon_dim=\"lon\",\n", - " lat_dim=\"lat\",\n", - ")\n", - "\n", - "qg_div_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" + "qg_div_data_500_202201 = ecl.calc_divergence(\n u_data=geostrophic_wind_data_500_202201.ug,\n v_data=geostrophic_wind_data_500_202201.vg,\n lon_dim=\"lon\",\n lat_dim=\"lat\",\n)\n\nqg_div_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Actual divergence\n", - "\n" + "Actual divergence\n\n" ] }, { @@ -524,25 +339,14 @@ }, "outputs": [], "source": [ - "div_data_500_202201 = ecl.calc_divergence(\n", - " u_data=uvdata_500_202201[\"uwnd\"],\n", - " v_data=uvdata_500_202201[\"vwnd\"],\n", - " lon_dim=\"lon\",\n", - " lat_dim=\"lat\",\n", - ")\n", - "\n", - "div_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" + "div_data_500_202201 = ecl.calc_divergence(\n u_data=uvdata_500_202201[\"uwnd\"],\n v_data=uvdata_500_202201[\"vwnd\"],\n lon_dim=\"lon\",\n lat_dim=\"lat\",\n)\n\ndiv_data_500_202201.sel(lat=slice(20, 80)).plot.contourf(levels=21)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Of course, in addition to the built-in finite difference method, the spherical harmonic function mothod can be solved, but you must ensure that it is **Global** and **Regular or Gaussian grid** type data.\n", - "\n", - "- :py:func:`easyclimate.windspharm.calc_relative_vorticity `: calculate the relative vorticity term with the spherical harmonic function mothod.\n", - "- :py:func:`easyclimate.windspharm.calc_divergence `: calculate the horizontal divergence term with the spherical harmonic function mothod.\n", - "\n" + "Of course, in addition to the built-in finite difference method, the spherical harmonic function mothod can be solved, but you must ensure that it is **Global** and **Regular or Gaussian grid** type data.\n\n- :py:func:`easyclimate.windspharm.calc_relative_vorticity `: calculate the relative vorticity term with the spherical harmonic function mothod.\n- :py:func:`easyclimate.windspharm.calc_divergence `: calculate the horizontal divergence term with the spherical harmonic function mothod.\n\n" ] }, { @@ -553,14 +357,7 @@ }, "outputs": [], "source": [ - "vor_data_500_202201_windspharm = ecl.windspharm.calc_relative_vorticity(\n", - " u_data=uvdata_500_202201[\"uwnd\"],\n", - " v_data=uvdata_500_202201[\"vwnd\"],\n", - ")\n", - "\n", - "vor_data_500_202201_windspharm.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(\n", - " levels=21\n", - ")" + "vor_data_500_202201_windspharm = ecl.windspharm.calc_relative_vorticity(\n u_data=uvdata_500_202201[\"uwnd\"],\n v_data=uvdata_500_202201[\"vwnd\"],\n)\n\nvor_data_500_202201_windspharm.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(\n levels=21\n)" ] }, { @@ -571,36 +368,21 @@ }, "outputs": [], "source": [ - "div_data_500_202201_windspharm = ecl.windspharm.calc_divergence(\n", - " u_data=uvdata_500_202201[\"uwnd\"],\n", - " v_data=uvdata_500_202201[\"vwnd\"],\n", - ")\n", - "\n", - "div_data_500_202201_windspharm.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(\n", - " levels=21\n", - ")" + "div_data_500_202201_windspharm = ecl.windspharm.calc_divergence(\n u_data=uvdata_500_202201[\"uwnd\"],\n v_data=uvdata_500_202201[\"vwnd\"],\n)\n\ndiv_data_500_202201_windspharm.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(\n levels=21\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Generally speaking, the calculation results of the finite difference method and the spherical harmonic function method are similar. The former does not require global regional data, but the calculation results of the latter are more accurate for high latitude regions.\n", - "\n" + "Generally speaking, the calculation results of the finite difference method and the spherical harmonic function method are similar. The former does not require global regional data, but the calculation results of the latter are more accurate for high latitude regions.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Advection\n", - "[Advection](https://glossary.ametsoc.org/wiki/Advection)_ is the process of transport of an atmospheric property solely by the mass motion (velocity field) of the atmosphere; also, the rate of change of the value of the advected property at a given point.\n", - "\n", - "For zonal advection, we can calculate as follows.\n", - "\n", - "\\begin{align}-u \\frac{\\partial T}{\\partial x}\\end{align}\n", - "\n", - "\n" + "## Advection\n[Advection](https://glossary.ametsoc.org/wiki/Advection)_ is the process of transport of an atmospheric property solely by the mass motion (velocity field) of the atmosphere; also, the rate of change of the value of the advected property at a given point.\n\nFor zonal advection, we can calculate as follows.\n\n\\begin{align}-u \\frac{\\partial T}{\\partial x}\\end{align}\n\n\n" ] }, { @@ -611,22 +393,14 @@ }, "outputs": [], "source": [ - "u_advection_500_202201 = ecl.calc_u_advection(\n", - " u_data=uvdata_500_202201[\"uwnd\"], temper_data=temp_data_500_202201\n", - ")\n", - "\n", - "u_advection_500_202201.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(levels=21)" + "u_advection_500_202201 = ecl.calc_u_advection(\n u_data=uvdata_500_202201[\"uwnd\"], temper_data=temp_data_500_202201\n)\n\nu_advection_500_202201.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(levels=21)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Similarly, the meridional advection can acquire as follows.\n", - "\n", - "\\begin{align}-v \\frac{\\partial T}{\\partial y}\\end{align}\n", - "\n", - "\n" + "Similarly, the meridional advection can acquire as follows.\n\n\\begin{align}-v \\frac{\\partial T}{\\partial y}\\end{align}\n\n\n" ] }, { @@ -637,31 +411,14 @@ }, "outputs": [], "source": [ - "v_advection_500_202201 = ecl.calc_v_advection(\n", - " v_data=uvdata_500_202201[\"vwnd\"], temper_data=temp_data_500_202201\n", - ")\n", - "\n", - "v_advection_500_202201.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(levels=21)" + "v_advection_500_202201 = ecl.calc_v_advection(\n v_data=uvdata_500_202201[\"vwnd\"], temper_data=temp_data_500_202201\n)\n\nv_advection_500_202201.sortby(\"lat\").sel(lat=slice(20, 80)).plot.contourf(levels=21)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Water Flux\n", - "\n", - "- :py:func:`easyclimate.calc_horizontal_water_flux `: calculate horizontal water vapor flux at each vertical level.\n", - "\n", - "\\begin{align}\\frac{1}{g} q \\mathbf{V} = \\frac{1}{g} (u q\\ \\mathbf{i} + vq\\ \\mathbf{j})\\end{align}\n", - "\n", - "- :py:func:`easyclimate.calc_vertical_water_flux `: calculate vertical water vapor flux.\n", - "\n", - "\\begin{align}-\\omega \\frac{q}{g}\\end{align}\n", - "\n", - "- :py:func:`easyclimate.calc_water_flux_top2surface_integral `: calculate the water vapor flux across the vertical level.\n", - "\n", - ":py:func:`easyclimate.calc_horizontal_water_flux ` can calculate the horizontal water flux of single layers.\n", - "\n" + "## Water Flux\n\n- :py:func:`easyclimate.calc_horizontal_water_flux `: calculate horizontal water vapor flux at each vertical level.\n\n\\begin{align}\\frac{1}{g} q \\mathbf{V} = \\frac{1}{g} (u q\\ \\mathbf{i} + vq\\ \\mathbf{j})\\end{align}\n\n- :py:func:`easyclimate.calc_vertical_water_flux `: calculate vertical water vapor flux.\n\n\\begin{align}-\\omega \\frac{q}{g}\\end{align}\n\n- :py:func:`easyclimate.calc_water_flux_top2surface_integral `: calculate the water vapor flux across the vertical level.\n\n:py:func:`easyclimate.calc_horizontal_water_flux ` can calculate the horizontal water flux of single layers.\n\n" ] }, { @@ -672,19 +429,14 @@ }, "outputs": [], "source": [ - "ecl.calc_horizontal_water_flux(\n", - " specific_humidity_data=q_data,\n", - " u_data=uvdata.uwnd,\n", - " v_data=uvdata.vwnd,\n", - ")" + "ecl.calc_horizontal_water_flux(\n specific_humidity_data=q_data,\n u_data=uvdata.uwnd,\n v_data=uvdata.vwnd,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The whole layer integral needs to consider the function :py:func:`easyclimate.calc_water_flux_top2surface_integral ` to calculate.\n", - "\n" + "The whole layer integral needs to consider the function :py:func:`easyclimate.calc_water_flux_top2surface_integral ` to calculate.\n\n" ] }, { @@ -695,25 +447,14 @@ }, "outputs": [], "source": [ - "water_flux_top2surface_integral = ecl.calc_water_flux_top2surface_integral(\n", - " specific_humidity_data=q_data,\n", - " u_data=u_data,\n", - " v_data=v_data,\n", - " surface_pressure_data=msl_data,\n", - " surface_pressure_data_units=\"millibars\",\n", - " vertical_dim=\"level\",\n", - " vertical_dim_units=\"hPa\",\n", - ")\n", - "\n", - "water_flux_top2surface_integral" + "water_flux_top2surface_integral = ecl.calc_water_flux_top2surface_integral(\n specific_humidity_data=q_data,\n u_data=u_data,\n v_data=v_data,\n surface_pressure_data=msl_data,\n surface_pressure_data_units=\"millibars\",\n vertical_dim=\"level\",\n vertical_dim_units=\"hPa\",\n)\n\nwater_flux_top2surface_integral" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Extracting the entire layer water vapor flux at mid and low latitudes at the 0th time level.\n", - "\n" + "Extracting the entire layer water vapor flux at mid and low latitudes at the 0th time level.\n\n" ] }, { @@ -724,12 +465,7 @@ }, "outputs": [], "source": [ - "draw_water_flux = (\n", - " water_flux_top2surface_integral.isel(time=0)\n", - " .thin(lon=3, lat=3)\n", - " .sel(lat=slice(-60, 60))\n", - ")\n", - "draw_pr = pr_data.isel(time=0).sel(lat=slice(-60, 60))" + "draw_water_flux = (\n water_flux_top2surface_integral.isel(time=0)\n .thin(lon=3, lat=3)\n .sel(lat=slice(-60, 60))\n)\ndraw_pr = pr_data.isel(time=0).sel(lat=slice(-60, 60))" ] }, { @@ -740,48 +476,14 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "draw_water_flux.plot.quiver(\n", - " ax=ax,\n", - " u=\"qu\",\n", - " v=\"qv\",\n", - " x=\"lon\",\n", - " y=\"lat\",\n", - " transform=ccrs.PlateCarree(),\n", - " zorder=2,\n", - ")\n", - "\n", - "draw_pr.plot.contourf(\n", - " ax=ax,\n", - " transform=ccrs.PlateCarree(),\n", - " levels=21,\n", - " cmap=\"Greens\",\n", - " zorder=1,\n", - " cbar_kwargs={\"location\": \"bottom\"},\n", - " vmax=20,\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\ndraw_water_flux.plot.quiver(\n ax=ax,\n u=\"qu\",\n v=\"qv\",\n x=\"lon\",\n y=\"lat\",\n transform=ccrs.PlateCarree(),\n zorder=2,\n)\n\ndraw_pr.plot.contourf(\n ax=ax,\n transform=ccrs.PlateCarree(),\n levels=21,\n cmap=\"Greens\",\n zorder=1,\n cbar_kwargs={\"location\": \"bottom\"},\n vmax=20,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Water Vapor Flux Divergence\n", - "\n", - "Water vapor flux divergence represents the convergence and divergence of water vapor. There are also two built-in functions to calculate the results of single-layers and whole-layer integration respectively.\n", - "\n", - "- :py:func:`easyclimate.calc_divergence_watervaporflux `: calculate water vapor flux divergence at each vertical level.\n", - "\n", - "\\begin{align}\\nabla \\left( \\frac{1}{g} q \\mathbf{V} \\right) = \\frac{1}{g} \\nabla \\cdot \\left( q \\mathbf{V} \\right)\\end{align}\n", - "\n", - "- :py:func:`easyclimate.calc_divergence_watervaporflux_top2surface_integral `: calculate water vapor flux divergence across the vertical level.\n", - "\n" + "## Water Vapor Flux Divergence\n\nWater vapor flux divergence represents the convergence and divergence of water vapor. There are also two built-in functions to calculate the results of single-layers and whole-layer integration respectively.\n\n- :py:func:`easyclimate.calc_divergence_watervaporflux `: calculate water vapor flux divergence at each vertical level.\n\n\\begin{align}\\nabla \\left( \\frac{1}{g} q \\mathbf{V} \\right) = \\frac{1}{g} \\nabla \\cdot \\left( q \\mathbf{V} \\right)\\end{align}\n\n- :py:func:`easyclimate.calc_divergence_watervaporflux_top2surface_integral `: calculate water vapor flux divergence across the vertical level.\n\n" ] }, { @@ -792,28 +494,14 @@ }, "outputs": [], "source": [ - "divergence_watervaporflux_top2surface_integral = (\n", - " ecl.calc_divergence_watervaporflux_top2surface_integral(\n", - " specific_humidity_data=q_data,\n", - " u_data=u_data,\n", - " v_data=v_data,\n", - " surface_pressure_data=msl_data,\n", - " surface_pressure_data_units=\"millibars\",\n", - " specific_humidity_data_units=\"grams/kg\",\n", - " vertical_dim=\"level\",\n", - " vertical_dim_units=\"hPa\",\n", - " )\n", - ")\n", - "\n", - "divergence_watervaporflux_top2surface_integral" + "divergence_watervaporflux_top2surface_integral = (\n ecl.calc_divergence_watervaporflux_top2surface_integral(\n specific_humidity_data=q_data,\n u_data=u_data,\n v_data=v_data,\n surface_pressure_data=msl_data,\n surface_pressure_data_units=\"millibars\",\n specific_humidity_data_units=\"grams/kg\",\n vertical_dim=\"level\",\n vertical_dim_units=\"hPa\",\n )\n)\n\ndivergence_watervaporflux_top2surface_integral" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Extracting the entire layer water vapor flux at mid and low latitudes at the 0th time level.\n", - "\n" + "Extracting the entire layer water vapor flux at mid and low latitudes at the 0th time level.\n\n" ] }, { @@ -824,9 +512,7 @@ }, "outputs": [], "source": [ - "draw_data = divergence_watervaporflux_top2surface_integral.isel(time=0).sel(\n", - " lat=slice(-60, 60)\n", - ")" + "draw_data = divergence_watervaporflux_top2surface_integral.isel(time=0).sel(\n lat=slice(-60, 60)\n)" ] }, { @@ -837,16 +523,7 @@ }, "outputs": [], "source": [ - "fig, ax = plt.subplots(\n", - " subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n", - ")\n", - "\n", - "ax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\n", - "ax.coastlines(edgecolor=\"black\", linewidths=0.5)\n", - "\n", - "draw_data.plot.contourf(\n", - " ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n", - ")" + "fig, ax = plt.subplots(\n subplot_kw={\"projection\": ccrs.PlateCarree(central_longitude=180)}\n)\n\nax.gridlines(draw_labels=[\"bottom\", \"left\"], color=\"grey\", alpha=0.5, linestyle=\"--\")\nax.coastlines(edgecolor=\"black\", linewidths=0.5)\n\ndraw_data.plot.contourf(\n ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={\"location\": \"bottom\"}, levels=21\n)" ] } ], @@ -866,7 +543,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/docs/example/plot_interp.ipynb b/docs/example/plot_interp.ipynb index 4cae2c0c..ae7bd4f5 100644 --- a/docs/example/plot_interp.ipynb +++ b/docs/example/plot_interp.ipynb @@ -143,6 +143,114 @@ "source": [ "fig, ax = plt.subplots(1, 2, figsize=(12, 5))\n\nu_data.sel(level=500).isel(time=0).plot(ax=ax[0])\nax[0].set_title(\"Before\", size=20)\n\nregriding_data.sel(level=500).isel(time=0).plot(ax=ax[1])\nax[1].set_title(\"After\", size=20)" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interpolation from model layers to altitude layers\nSuppose the following data are available\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "uwnd_data = xr.DataArray(\n np.array([-15.080393 , -10.749071 , -4.7920494, -2.3813725, -1.4152431,\n -0.6465206, -7.8181705, -14.718096 , -14.65539 , -14.948015 ,\n -13.705519 , -11.443476 , -8.865583 , -7.9528713, -8.329103 ,\n -7.2445316, -6.7150173, -5.5189686, -4.139448 , -3.2731838,\n -2.2931194, -1.0041752, -1.8983078, -2.3674374, -2.8203106,\n -3.2940865, -3.526329 , -3.8654022, -4.164995 , -4.2834396,\n -4.2950516, -4.3438225, -4.3716908, -4.7688255, -4.6155453,\n -4.5528393, -4.4831676, -4.385626 , -4.2950516, -4.0953217]),\n dims = ('model_level',),\n coords = {'model_level': np.array([ 36, 44, 51, 56, 60, 63, 67, 70, 73, 75, 78, 81, 83, 85,\n 88, 90, 92, 94, 96, 98, 100, 102, 104, 105, 107, 109, 111, 113,\n 115, 117, 119, 122, 124, 126, 129, 131, 133, 135, 136, 137])\n }\n)\n\np_data = xr.DataArray(\n np.array([ 2081.4756, 3917.6995, 6162.6455, 8171.3506, 10112.652 ,\n 11811.783 , 14447.391 , 16734.607 , 19317.787 , 21218.21 ,\n 24357.875 , 27871.277 , 30436.492 , 33191.027 , 37698.96 ,\n 40969.438 , 44463.73 , 48191.92 , 52151.29 , 56291.098 ,\n 60525.63 , 64770.367 , 68943.69 , 70979.66 , 74908.17 ,\n 78599.67 , 82012.95 , 85122.69 , 87918.29 , 90401.77 ,\n 92584.94 , 95338.72 , 96862.08 , 98165.305 , 99763.38 ,\n 100626.21 , 101352.69 , 101962.28 , 102228.875 , 102483.055 ]),\n dims = ('model_level',),\n coords = {'model_level': np.array([ 36, 44, 51, 56, 60, 63, 67, 70, 73, 75, 78, 81, 83, 85,\n 88, 90, 92, 94, 96, 98, 100, 102, 104, 105, 107, 109, 111, 113,\n 115, 117, 119, 122, 124, 126, 129, 131, 133, 135, 136, 137])\n }\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we interpolate the data located on the mode plane to the isobaric plane. \nNote that the units of the given isobaric surface are consistent with `pressure_data`.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "result = ecl.interp.interp1d_vertical_model2pressure(\n pressure_data = p_data,\n variable_data = uwnd_data,\n vertical_input_dim = 'model_level',\n vertical_output_dim = 'plev',\n vertical_output_level = np.array([100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 20000, 10000])\n)\nresult" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Simply calibrate the interpolation effect.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\nax.plot(p_data, uwnd_data, label = 'Original')\nax.plot(result.plev, result.data, 'o', label = 'Interpolated')\nax.invert_xaxis()\nax.set_xlabel('Pressure [Pa]')\nax.set_ylabel('Zonal Wind [m/s]')\nplt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interpolation from pressure layers to altitude layers\nSuppose the following data are available\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "z_data = xr.DataArray(\n np.array([ 214.19354, 841.6774 , 1516.871 , 3055.7097 , 4260.5806 ,\n 5651.4194 , 7288.032 , 9288.193 , 10501.097 , 11946.71 ,\n 13762.322 , 16233.451 , 18370.902 , 20415.227 , 23619.033 ,\n 26214.322 , 30731.807 ]),\n dims=('level'),\n coords={'level': np.array([1000., 925., 850., 700., 600., 500., 400., 300., 250., 200.,\n 150., 100., 70., 50., 30., 20., 10.])}\n)\nuwnd_data = xr.DataArray(\n np.array([-2.3200073, -3.5700073, -2.5800018, 8.080002 , 14.059998 ,\n 22.119995 , 33.819992 , 49.339996 , 57.86 , 64.009995 ,\n 62.940002 , 49.809998 , 31.160004 , 16.59999 , 10.300003 ,\n 10.459991 , 9.880005 ]),\n dims=('level'), \n coords={'level': np.array([1000., 925., 850., 700., 600., 500., 400., 300., 250., 200.,\n 150., 100., 70., 50., 30., 20., 10.])} \n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we need to interpolate the zonal wind data (located on the isobaric surface) to the altitude layers.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "target_heights = np.linspace(0, 10000, 101)\n\nresult = ecl.interp.interp1d_vertical_pressure2altitude(\n z_data = z_data,\n variable_data = uwnd_data,\n target_heights = target_heights,\n vertical_input_dim = 'level',\n vertical_output_dim = 'altitude',\n)\nresult" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can check the interpolation results.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plt.plot(z_data[:9], uwnd_data[:9], 'o', label = 'Original')\nplt.plot(result.altitude, result.data, label = 'Interpolated')\nplt.xlabel('Altitude [m]')\nplt.ylabel('Zonal Wind [m/s]')\nplt.legend()" + ] } ], "metadata": { @@ -161,7 +269,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/docs/example/plot_taylor_diagram.ipynb b/docs/example/plot_taylor_diagram.ipynb index 80004fe4..e582df9b 100644 --- a/docs/example/plot_taylor_diagram.ipynb +++ b/docs/example/plot_taylor_diagram.ipynb @@ -197,7 +197,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/docs/example/plot_time_scale_average.ipynb b/docs/example/plot_time_scale_average.ipynb index af77c421..8c2fce56 100644 --- a/docs/example/plot_time_scale_average.ipynb +++ b/docs/example/plot_time_scale_average.ipynb @@ -143,7 +143,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.15" } }, "nbformat": 4, diff --git a/docs/source/dynamic_docs/plot_interp.py b/docs/source/dynamic_docs/plot_interp.py index fb45c576..f9b8825a 100644 --- a/docs/source/dynamic_docs/plot_interp.py +++ b/docs/source/dynamic_docs/plot_interp.py @@ -90,3 +90,344 @@ regriding_data.sel(level=500).isel(time=0).plot(ax=ax[1]) ax[1].set_title("After", size=20) + +# %% +# Interpolation from model layers to altitude layers +# ------------------------------------ +# Suppose the following data are available + +uwnd_data = xr.DataArray( + np.array( + [ + -15.080393, + -10.749071, + -4.7920494, + -2.3813725, + -1.4152431, + -0.6465206, + -7.8181705, + -14.718096, + -14.65539, + -14.948015, + -13.705519, + -11.443476, + -8.865583, + -7.9528713, + -8.329103, + -7.2445316, + -6.7150173, + -5.5189686, + -4.139448, + -3.2731838, + -2.2931194, + -1.0041752, + -1.8983078, + -2.3674374, + -2.8203106, + -3.2940865, + -3.526329, + -3.8654022, + -4.164995, + -4.2834396, + -4.2950516, + -4.3438225, + -4.3716908, + -4.7688255, + -4.6155453, + -4.5528393, + -4.4831676, + -4.385626, + -4.2950516, + -4.0953217, + ] + ), + dims=("model_level",), + coords={ + "model_level": np.array( + [ + 36, + 44, + 51, + 56, + 60, + 63, + 67, + 70, + 73, + 75, + 78, + 81, + 83, + 85, + 88, + 90, + 92, + 94, + 96, + 98, + 100, + 102, + 104, + 105, + 107, + 109, + 111, + 113, + 115, + 117, + 119, + 122, + 124, + 126, + 129, + 131, + 133, + 135, + 136, + 137, + ] + ) + }, +) + +p_data = xr.DataArray( + np.array( + [ + 2081.4756, + 3917.6995, + 6162.6455, + 8171.3506, + 10112.652, + 11811.783, + 14447.391, + 16734.607, + 19317.787, + 21218.21, + 24357.875, + 27871.277, + 30436.492, + 33191.027, + 37698.96, + 40969.438, + 44463.73, + 48191.92, + 52151.29, + 56291.098, + 60525.63, + 64770.367, + 68943.69, + 70979.66, + 74908.17, + 78599.67, + 82012.95, + 85122.69, + 87918.29, + 90401.77, + 92584.94, + 95338.72, + 96862.08, + 98165.305, + 99763.38, + 100626.21, + 101352.69, + 101962.28, + 102228.875, + 102483.055, + ] + ), + dims=("model_level",), + coords={ + "model_level": np.array( + [ + 36, + 44, + 51, + 56, + 60, + 63, + 67, + 70, + 73, + 75, + 78, + 81, + 83, + 85, + 88, + 90, + 92, + 94, + 96, + 98, + 100, + 102, + 104, + 105, + 107, + 109, + 111, + 113, + 115, + 117, + 119, + 122, + 124, + 126, + 129, + 131, + 133, + 135, + 136, + 137, + ] + ) + }, +) + +# %% +# Now we interpolate the data located on the mode plane to the isobaric plane. +# Note that the units of the given isobaric surface are consistent with `pressure_data`. + +result = ecl.interp.interp1d_vertical_model2pressure( + pressure_data=p_data, + variable_data=uwnd_data, + vertical_input_dim="model_level", + vertical_output_dim="plev", + vertical_output_level=np.array( + [100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 20000, 10000] + ), +) +result + +# %% +# Simply calibrate the interpolation effect. + +fig, ax = plt.subplots() +ax.plot(p_data, uwnd_data, label="Original") +ax.plot(result.plev, result.data, "o", label="Interpolated") +ax.invert_xaxis() +ax.set_xlabel("Pressure [Pa]") +ax.set_ylabel("Zonal Wind [m/s]") +plt.legend() + +# %% +# Interpolation from pressure layers to altitude layers +# ------------------------------------ +# Suppose the following data are available + +z_data = xr.DataArray( + np.array( + [ + 214.19354, + 841.6774, + 1516.871, + 3055.7097, + 4260.5806, + 5651.4194, + 7288.032, + 9288.193, + 10501.097, + 11946.71, + 13762.322, + 16233.451, + 18370.902, + 20415.227, + 23619.033, + 26214.322, + 30731.807, + ] + ), + dims=("level"), + coords={ + "level": np.array( + [ + 1000.0, + 925.0, + 850.0, + 700.0, + 600.0, + 500.0, + 400.0, + 300.0, + 250.0, + 200.0, + 150.0, + 100.0, + 70.0, + 50.0, + 30.0, + 20.0, + 10.0, + ] + ) + }, +) +uwnd_data = xr.DataArray( + np.array( + [ + -2.3200073, + -3.5700073, + -2.5800018, + 8.080002, + 14.059998, + 22.119995, + 33.819992, + 49.339996, + 57.86, + 64.009995, + 62.940002, + 49.809998, + 31.160004, + 16.59999, + 10.300003, + 10.459991, + 9.880005, + ] + ), + dims=("level"), + coords={ + "level": np.array( + [ + 1000.0, + 925.0, + 850.0, + 700.0, + 600.0, + 500.0, + 400.0, + 300.0, + 250.0, + 200.0, + 150.0, + 100.0, + 70.0, + 50.0, + 30.0, + 20.0, + 10.0, + ] + ) + }, +) + +# %% +# Then we need to interpolate the zonal wind data (located on the isobaric surface) to the altitude layers. + +target_heights = np.linspace(0, 10000, 101) + +result = ecl.interp.interp1d_vertical_pressure2altitude( + z_data=z_data, + variable_data=uwnd_data, + target_heights=target_heights, + vertical_input_dim="level", + vertical_output_dim="altitude", +) +result + +# %% +# Now we can check the interpolation results. +plt.plot(z_data[:9], uwnd_data[:9], "o", label="Original") +plt.plot(result.altitude, result.data, label="Interpolated") +plt.xlabel("Altitude [m]") +plt.ylabel("Zonal Wind [m/s]") +plt.legend() diff --git a/src/easyclimate/interp/__init__.py b/src/easyclimate/interp/__init__.py index 6dbd21c0..7abbaf87 100644 --- a/src/easyclimate/interp/__init__.py +++ b/src/easyclimate/interp/__init__.py @@ -1,3 +1,4 @@ from .mesh2mesh import * from .point2mesh import * from .modellevel2pressure import * +from .pressure2altitude import * diff --git a/src/easyclimate/interp/pressure2altitude.py b/src/easyclimate/interp/pressure2altitude.py new file mode 100644 index 00000000..520b1d23 --- /dev/null +++ b/src/easyclimate/interp/pressure2altitude.py @@ -0,0 +1,97 @@ +""" +Interpolate from pressure layer to altitude layers +""" + +from __future__ import annotations + +import numpy as np +import xarray as xr + +__all__ = ["interp1d_vertical_pressure2altitude"] + + +def interp1d_vertical_pressure2altitude( + z_data: xr.DataArray, + variable_data: xr.DataArray, + target_heights: np.array, + vertical_input_dim: str, + vertical_output_dim: str, + kind: str = "cubic", + bounds_error=None, + fill_value="extrapolate", + assume_sorted: bool = False, +) -> xr.DataArray: + """ + Interpolating variables from the pressure levels (e.g., sigma vertical coordinate) to the altitude levels by 1-D function. + + z_data: :py:class:`xarray.DataArray`. + The vertical atmospheric geopotential height on each pressure level. + variable_data: :py:class:`xarray.DataArray`. + variable to be interpolated on model levels. + + .. note:: + The shape of `z_data` `z_data` and `variable_data` shuold be the same. + + target_heights: :py:class:`numpy.array`. + Altitude interpolation sampling range. e.g., `np.linspace(0, 15000, 100)`. + + .. note:: + The unit of `target_heights` shuold be same as `z_data`, e.g., the unit of `target_heights` is `meter`, and the unit of `z_data` shuold be `meter`. + + vertical_input_dim: :py:class:`str `. + The name of the standard pressure levels dimension, often assigned the value `'plev'` (Pa) or `'lev'` (hPa). + + .. note:: + The `vertical_input_dim` `z_data` and `variable_data` shuold be the same. + + vertical_output_dim: :py:class:`str `. + The name of the altitude dimension, often assigned the value `'altitude'`. + kind: :py:class:`str ` or :py:class:`int `, default: `'cubic'`. + Specifies the kind of interpolation as a string or as an integer specifying the order of the spline interpolator to use. + The string has to be one of 'linear', 'nearest', 'nearest-up', 'zero', 'slinear', 'quadratic', 'cubic', 'previous', + or 'next'. 'zero', 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation of zeroth, first, second or third order; + 'previous' and 'next' simply return the previous or next value of the point; + 'nearest-up' and 'nearest' differ when interpolating half-integers (e.g. 0.5, 1.5) in that 'nearest-up' rounds up and 'nearest' rounds down. + bounds_error: :py:class:`bool `, default: `None`. + If True, a ValueError is raised any time interpolation is attempted on a value outside of the range of `pressure_data` (where extrapolation is necessary). + If False, out of bounds values are assigned fill_value. By default, an error is raised unless `fill_value="extrapolate"`. + fill_value: array-like or (array-like, array_like) or "extrapolate", default: `extrapolate`. + - If a ndarray (or float), this value will be used to fill in for requested points outside of the data range. If not provided, then the default is NaN. The array-like must broadcast properly to the dimensions of the non-interpolation axes. + - If a two-element tuple, then the first element is used as a fill value for x_new < x[0] and the second element is used for x_new > x[-1]. Anything that is not a 2-element tuple (e.g., list or ndarray, regardless of shape) is taken to be a single array-like argument meant to be used for both bounds as below, above = fill_value, fill_value. Using a two-element tuple or ndarray requires bounds_error=False. + - If "extrapolate", then points outside the data range will be extrapolated. + assume_sortedbool: :py:class:`bool `, default: `False`. + If `False`, values of x can be in any order and they are sorted first. If `True`, `pressure_data` has to be an array of monotonically increasing values. + + .. seealso:: + - `scipy.interpolate.interp1d `__ + """ + from scipy import interpolate + + def _interp_height(z_data, variable_data, target_heights=target_heights): + object_interp = interpolate.interp1d( + z_data, + variable_data, + kind=kind, + bounds_error=bounds_error, + fill_value=fill_value, + assume_sorted=assume_sorted, + ) + result_interp = object_interp(target_heights) + return result_interp + + interped_data = xr.apply_ufunc( + _interp_height, + z_data, + variable_data, + input_core_dims=[[vertical_input_dim], [vertical_input_dim]], + output_core_dims=[[vertical_output_dim]], + output_dtypes=["float64"], + dask="parallelized", + vectorize=True, + dask_gufunc_kwargs={ + "output_sizes": {vertical_output_dim: len(target_heights)}, + "allow_rechunk": True, + }, + ) + interped_data = interped_data.assign_coords({vertical_output_dim: target_heights}) + return interped_data diff --git a/test/test_interp_pressure2altitude.py b/test/test_interp_pressure2altitude.py new file mode 100644 index 00000000..b99db35e --- /dev/null +++ b/test/test_interp_pressure2altitude.py @@ -0,0 +1,221 @@ +""" +pytest for interp.pressure2altitude.py +""" + +import pytest + +import easyclimate as ecl +import numpy as np +import xarray as xr + +z_data = xr.DataArray( + np.array( + [ + 214.19354, + 841.6774, + 1516.871, + 3055.7097, + 4260.5806, + 5651.4194, + 7288.032, + 9288.193, + 10501.097, + 11946.71, + 13762.322, + 16233.451, + 18370.902, + 20415.227, + 23619.033, + 26214.322, + 30731.807, + ] + ), + dims=("level"), + coords={ + "level": np.array( + [ + 1000.0, + 925.0, + 850.0, + 700.0, + 600.0, + 500.0, + 400.0, + 300.0, + 250.0, + 200.0, + 150.0, + 100.0, + 70.0, + 50.0, + 30.0, + 20.0, + 10.0, + ] + ) + }, +) +uwnd_data = xr.DataArray( + np.array( + [ + -2.3200073, + -3.5700073, + -2.5800018, + 8.080002, + 14.059998, + 22.119995, + 33.819992, + 49.339996, + 57.86, + 64.009995, + 62.940002, + 49.809998, + 31.160004, + 16.59999, + 10.300003, + 10.459991, + 9.880005, + ] + ), + dims=("level"), + coords={ + "level": np.array( + [ + 1000.0, + 925.0, + 850.0, + 700.0, + 600.0, + 500.0, + 400.0, + 300.0, + 250.0, + 200.0, + 150.0, + 100.0, + 70.0, + 50.0, + 30.0, + 20.0, + 10.0, + ] + ) + }, +) + + +def test_interp1d_vertical_pressure2altitude(): + target_heights = np.linspace(0, 10000, 101) + result_data = ecl.interp.interp1d_vertical_pressure2altitude( + z_data=z_data, + variable_data=uwnd_data, + target_heights=target_heights, + vertical_input_dim="level", + vertical_output_dim="altitude", + ) + refer_data = np.array( + [ + -1.62144403, + -1.95839842, + -2.27662081, + -2.5715752, + -2.83872558, + -3.07353596, + -3.27147036, + -3.42799276, + -3.53856718, + -3.59865762, + -3.60372809, + -3.54924259, + -3.43066512, + -3.24345969, + -2.9830903, + -2.64502095, + -2.22597055, + -1.73105429, + -1.16883255, + -0.54787619, + 0.12324392, + 0.83595692, + 1.58169195, + 2.35187814, + 3.13794463, + 3.93132056, + 4.72343507, + 5.50571728, + 6.26959634, + 7.00650139, + 7.70786156, + 8.3653115, + 8.97676994, + 9.54745, + 10.08297379, + 10.58896343, + 11.07104104, + 11.53482871, + 11.98594858, + 12.43002274, + 12.87267333, + 13.31952244, + 13.77619219, + 14.24823926, + 14.73858676, + 15.24668401, + 15.77174282, + 16.31297497, + 16.86959227, + 17.4408065, + 18.02582947, + 18.62387297, + 19.23414879, + 19.85586874, + 20.4882446, + 21.13048817, + 21.78181125, + 22.44143667, + 23.10885889, + 23.78385428, + 24.46621225, + 25.15572225, + 25.85217371, + 26.55535608, + 27.26505878, + 27.98107125, + 28.70318293, + 29.43118325, + 30.16486166, + 30.90400758, + 31.64841045, + 32.39785972, + 33.15214481, + 33.91105509, + 34.6743274, + 35.44155116, + 36.21229011, + 36.98610801, + 37.76256861, + 38.54123568, + 39.32167296, + 40.1034442, + 40.88611317, + 41.66924361, + 42.45239929, + 43.23514394, + 44.01704134, + 44.79765523, + 45.57654936, + 46.3532875, + 47.12743339, + 47.8985508, + 48.66620346, + 49.42995489, + 50.18915079, + 50.94252301, + 51.68869596, + 52.42629411, + 53.15394187, + 53.8702637, + 54.57388404, + ] + ) + assert np.isclose(result_data, refer_data).all() From be417aa5a92dde2b833cc835f1bfe15d2d4a7c3e Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 20:45:05 +0800 Subject: [PATCH 2/8] fix: Fix whitespace and formatting issues in docs_readthedocs.yml --- .github/workflows/docs_readthedocs.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/docs_readthedocs.yml b/.github/workflows/docs_readthedocs.yml index d10da1df..c49734d7 100644 --- a/.github/workflows/docs_readthedocs.yml +++ b/.github/workflows/docs_readthedocs.yml @@ -61,7 +61,6 @@ jobs: with: path: ./docs/build/html - - name: Trigger RTDs build uses: dfm/rtds-action@v1 with: From fca5fec9a860b32786ceffac80dfccd4cab97a2e Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 20:52:59 +0800 Subject: [PATCH 3/8] feat: add docs seaborn rely --- docs/requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/requirements.txt b/docs/requirements.txt index ec2d3bf8..47dbf8dd 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -36,3 +36,4 @@ sphinx-gallery == 0.17.1 sphinx-autoapi == 3.3.2 sphinx-copybutton == 0.5.2 sphinx-book-theme == 1.1.3 +seaborn From 688dbbb379c4ccc49700fad23a2d1c5ddd870cec Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 21:29:16 +0800 Subject: [PATCH 4/8] Update: docs_readthedocs.yml --- .github/workflows/docs_readthedocs.yml | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/.github/workflows/docs_readthedocs.yml b/.github/workflows/docs_readthedocs.yml index c49734d7..2f94fe9c 100644 --- a/.github/workflows/docs_readthedocs.yml +++ b/.github/workflows/docs_readthedocs.yml @@ -52,18 +52,13 @@ jobs: uses: ad-m/github-push-action@v0.8.0 with: github_token: ${{ secrets.GITHUB_TOKEN }} - branch: gh-pages + branch: main-pages force: true directory: ./docs/build/html - - name: Upload artifact - uses: actions/upload-pages-artifact@v3 - with: - path: ./docs/build/html - - name: Trigger RTDs build - uses: dfm/rtds-action@v1 - with: - webhook_url: ${{ secrets.RTDS_WEBHOOK_URL }} - webhook_token: ${{ secrets.RTDS_WEBHOOK_TOKEN }} - commit_ref: ${{ github.ref }} + run: | + curl -X POST \ + -H "Content-Type: application/x-www-form-urlencoded" \ + -d "branches=dev&token=${{ secrets.RTDS_WEBHOOK_TOKEN }}" \ + ${{ secrets.RTDS_WEBHOOK_URL }} From a82695a48c918cd792cb84bb47fc125a5e4d3e46 Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 22:29:16 +0800 Subject: [PATCH 5/8] Update: docs_readthedocs.yml --- .github/workflows/docs_readthedocs.yml | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/.github/workflows/docs_readthedocs.yml b/.github/workflows/docs_readthedocs.yml index 2f94fe9c..f9d2e4bb 100644 --- a/.github/workflows/docs_readthedocs.yml +++ b/.github/workflows/docs_readthedocs.yml @@ -48,17 +48,14 @@ jobs: git config --local user.name "GitHub Action" git commit -m 'deploy' - - name: Force push to destination branch 🚀 - uses: ad-m/github-push-action@v0.8.0 + - uses: actions/upload-artifact@v3 with: - github_token: ${{ secrets.GITHUB_TOKEN }} - branch: main-pages - force: true - directory: ./docs/build/html + name: notebooks-for-${{ github.sha }} + path: docs/build/html - name: Trigger RTDs build - run: | - curl -X POST \ - -H "Content-Type: application/x-www-form-urlencoded" \ - -d "branches=dev&token=${{ secrets.RTDS_WEBHOOK_TOKEN }}" \ - ${{ secrets.RTDS_WEBHOOK_URL }} + uses: dfm/rtds-action@v1 + with: + webhook_url: ${{ secrets.RTDS_WEBHOOK_URL }} + webhook_token: ${{ secrets.RTDS_WEBHOOK_TOKEN }} + commit_ref: ${{ github.ref }} From 6f47d659e1cbf78ae858d422c1abe7b57423563c Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 22:29:49 +0800 Subject: [PATCH 6/8] Delete docs_readthedocs.yml --- .github/workflows/docs_readthedocs.yml | 61 -------------------------- 1 file changed, 61 deletions(-) delete mode 100644 .github/workflows/docs_readthedocs.yml diff --git a/.github/workflows/docs_readthedocs.yml b/.github/workflows/docs_readthedocs.yml deleted file mode 100644 index f9d2e4bb..00000000 --- a/.github/workflows/docs_readthedocs.yml +++ /dev/null @@ -1,61 +0,0 @@ -name: docs_readthedocs - -# execute this workflow automatically when a we push to master -on: - push: - branches: [ main, dev ] - -# Sets permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages -# https://github.com/ad-m/github-push-action/issues/52 -permissions: write-all - -jobs: - - build_docs_job: - runs-on: ubuntu-latest - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - - steps: - - name: Checkout 🛎️ - uses: actions/checkout@v4.2.2 - - - name: Set up Python - uses: actions/setup-python@v5.3.0 - with: - python-version: '3.10' - - - name: Install OptiPNG and Build 🔧 - run: | - sudo apt-get install -y optipng - - - name: Install dependencies 💻 - run: | - pip install -r docs/requirements.txt - - - name: make the sphinx docs - run: | - # make -C docs clean - make -C docs html - - - name: Init new repo in dist folder and commit generated files - run: | - cd docs/build/html/ - git init - touch .nojekyll - git add -A - git config --local user.email "action@github.com" - git config --local user.name "GitHub Action" - git commit -m 'deploy' - - - uses: actions/upload-artifact@v3 - with: - name: notebooks-for-${{ github.sha }} - path: docs/build/html - - - name: Trigger RTDs build - uses: dfm/rtds-action@v1 - with: - webhook_url: ${{ secrets.RTDS_WEBHOOK_URL }} - webhook_token: ${{ secrets.RTDS_WEBHOOK_TOKEN }} - commit_ref: ${{ github.ref }} From cef99c25a18cf017623926d4bec135c745fb8a1b Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 22:51:07 +0800 Subject: [PATCH 7/8] fix: rtd github download --- src/easyclimate/core/tutorial.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/easyclimate/core/tutorial.py b/src/easyclimate/core/tutorial.py index b2267133..46b0c6b1 100644 --- a/src/easyclimate/core/tutorial.py +++ b/src/easyclimate/core/tutorial.py @@ -157,6 +157,9 @@ def open_tutorial_dataset( logger = pooch.get_logger() logger.setLevel("WARNING") + # pass user-agent to allow rtd downloading from github + downloader = pooch.HTTPDownloader(headers={"User-Agent": "easyclimate"}) + cache_dir = _construct_cache_dir(cache_dir) if name in external_urls: url = external_urls[name] @@ -182,7 +185,11 @@ def open_tutorial_dataset( # retrieve the file filepath = pooch.retrieve( - url=url, known_hash=None, path=cache_dir, progressbar=progressbar + url=url, + known_hash=None, + path=cache_dir, + progressbar=progressbar, + downloader=downloader, ) if Path(filepath).suffix == ".nc" or Path(filepath).suffix == ".grib": From 6cc13cdb054b7f871d8366333692f314668ba63e Mon Sep 17 00:00:00 2001 From: Shenyulu <59901836+shenyulu@users.noreply.github.com> Date: Tue, 3 Dec 2024 23:03:24 +0800 Subject: [PATCH 8/8] fix: seaborn relay --- docs/requirements.txt | 2 +- release_requirements.txt | 1 + test_requirements.txt | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index 47dbf8dd..d1a2ae9e 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -26,6 +26,7 @@ metpy tcpypi netCDF4 h5netcdf +seaborn sphinx == 8.0.2 recommonmark == 0.7.1 @@ -36,4 +37,3 @@ sphinx-gallery == 0.17.1 sphinx-autoapi == 3.3.2 sphinx-copybutton == 0.5.2 sphinx-book-theme == 1.1.3 -seaborn diff --git a/release_requirements.txt b/release_requirements.txt index 7a0716e7..b887ee8e 100644 --- a/release_requirements.txt +++ b/release_requirements.txt @@ -26,3 +26,4 @@ metpy tcpypi netCDF4 h5netcdf +seaborn diff --git a/test_requirements.txt b/test_requirements.txt index 124702dc..618706d9 100644 --- a/test_requirements.txt +++ b/test_requirements.txt @@ -26,6 +26,7 @@ metpy tcpypi netCDF4 h5netcdf +seaborn pytest == 8.3.3 pytest-cov == 5.0.0