Pythonを使用して気候データのゾーン統計を実行する

May 08 2022
あなたの国の気候は他の国と比べてどのように変化していますか?気候変動は、私たち全員の未来を徐々に形作り始めています。今後数十年で、気温の上昇、干ばつが長引く、ハリケーンの強さなどが予想されます。

あなたの国の気候は他の国と比べてどのように変化していますか?

6月の気温はヨーロッパのほとんどの国で上昇しています(著者による画像)

気候変動は、私たち全員の未来を徐々に形作り始めています。今後数十年で、気温の上昇、干ばつが長引く、ハリケーンの強さなどが予想されます。

インターネット上でこのトピックに関する多くの議論を見つけることができるので、この記事では気候変動については議論しません。ただし、欧州委員会とコペルニクスからの無料データを調べて、気温が何年にもわたってどのように変化しているかを理解するのに役立てます。

この記事では、複数のデータセットに対してゾーン統計を実行します。月平均気温の異常を可視化し、複数の国の時系列を作成します。最後に、ヒートマップチャートを使用して調査結果を視覚化します。

TLDR:データの取得方法とセットアップ環境を気にしない場合は、[ゾーン統計の実行]セクションにジャンプして、記事の一部のコーディングを楽しんでください。

環境のセットアップ

Python環境をセットアップする

Anacondaを使用した新しいPython環境の作成から始めます。xarray、regionmask、geopandasなどの多くの地理空間ライブラリを使用します。また、cdsapiを使用してClimate Data Storeに接続し、データを取得します。

conda create -n climate_change python=3.8 pandas scipy cdsapi xarray geopandas cfgrib regionmask

conda activate climate_change

次に、ClimateDataStoreからデータをダウンロードするために使用するCDSPythonAPIを正しく設定します。私の最後の(そして最初の)記事に従って進んでください。

すべてのデータとipythonノートブックを配置する新しいフォルダーを作成してください。このフォルダーで以下のすべての手順を実行します。この記事の最後で共有するノートブックclimate_change.ipynbを作成しました。

使用するデータを取得する

以前の記事の知識を再利用し、ClimateDataStoreをもう一度探索します。

気候データを取得する

今日は、コペルニクスによる気候変動の評価のために、必須の気候変数を使用します。データセットは、地表気温、地表空気相対湿度降水量 などの複数の気候変数で構成されています。次に、異常(気候学から月平均を差し引いたもの)、気候学、または月平均の製品タイプを選択できます。

大量のデータがあるため、1979年から2022年までの44年間の1月と6月の気温を比較します(この記事は2022年5月に作成されているため、6月のデータセットは43個しかないことに注意してください)。最終的に87個のファイルが必要です。

PythonAPIリクエストは次のようになります。

import cdsapi
c = cdsapi.Client()
c.retrieve(
    'ecv-for-climate-change',
    {
        'format': 'zip',
        'variable': 'surface_air_temperature',
        'climate_reference_period': '1991_2020',
        'product_type': [
            'anomaly', 'monthly_mean',
        ],
        'time_aggregation': '1_month_mean',
        'year': [
            '1979', '1980', '1981',
            '1982', '1983', '1984',
            '1985', '1986', '1987',
            '1988', '1989', '1990',
            '1991', '1992', '1993',
            '1994', '1995', '1996',
            '1997', '1998', '1999',
            '2000', '2001', '2002',
            '2003', '2004', '2005',
            '2006', '2007', '2008',
            '2009', '2010', '2011',
            '2012', '2013', '2014',
            '2015', '2016', '2017',
            '2018', '2019', '2020',
            '2021', '2022',
        ],
        'month': ['01','06'],
        'origin': 'era5',
    },
    'download.zip')

Once the data is downloaded, you should unzip it to folder of your choice. In result, you should have 87 grib files placed in your desired folder.

Create/Get Shapefiles

To calculate zonal statistics, we want to calculate statistics over specific area. In our case, the specific areas are European countries. We will use NUTS shapefilea file containing geospatial information about countries of European Union.

We will use shapefile, that was already created by European Commission. However, if you are familiar with GIS systems, you can create one of your own.

Our shapefile dataset with description is available for download here. We will use this configuration and download it.

Configuration for downloading shapefile used in this article

Save it to preferred location, I strongly suggest it being the same directory as the one where we store python notebook.

The NUTS are a hierarchical system divided into 3 levels. NUTS 1: major socio-economic regions, NUTS 2: basic regions for the application of regional policies, NUTS 3: small regions for specific diagnoses.

We will use NUTS1 and in bonus example, we will work with NUTS2.

This is how our NUTS shapefile looks like in QGIS, adding some basic symbology to differentiate between countries:

NUTS dataset visualised in QGIS (Image by Author)

Perform Zonal Statistics

Since we have gathered all of the data that we need, we can continue with our second to last step — Performing Zonal Statistics. What we want to achieve is to compute mean temperature for every element (country) in our shapefile over the course of 44 years for January and 43 years for June, respectively. We will perform computation for product type Anomaly only (you can do Monthly mean on your own after reading this article). We will extract this information to Pandas DataFrame, which we will use to visualise our results in the end.

As a bonus, we will select one country and perform zonal statistics on different regions within the country.

Let’s dive into coding part of this article.

Explore Climate Dataset

import xarray as xr 
import glob
import pandas as pd
# Get list of files with anomaly data using glob
grib_files = glob.glob("data/ecv/temperature/*anomaly*.grib")
# Explore dataset
with xr.open_dataset(grib_files[0], engine="cfgrib") as data:
    ds = data
ds

      
                

ただし、これはxarray open_mfdataset()関数を使用して簡単に実行できます。この関数では、開くファイルのリストを指定したり、ネストと連結の次元を指定したりできます。この例では、 valid_timeでデータを連結します。これは、上記のデータセットの座標の下に表示されます。open_mfdataset()に適切な引数を設定すると、必要なすべてのファイルを3次元データセットに直接開くことができます。また、すべてのオープニングは、複数のプロセッサで並行して実行できます。

# Open multifile dataset with xarray, concat them on "time" 
# dimension
with xr.open_mfdataset(grib_files, combine="nested",concat_dim = "valid_time", parallel=True) as data:
    multi_ds = data
multi_ds

      
                

import geopandas as gpd
eu_countries_shp = "data/shapefiles/NUTS_RG_01M_2021_4326.shp"
eu_countries_gdf = gpd.read_file(eu_countries_shp)
eu_countries_gdf.head()

      
                

eu_countries_gdf = eu_countries_gdf[eu_countries_gdf["LEVL_CODE"] == 0]

リージョン全体のゾーン平均を計算するには、最初にライブラリregionmaskを使用して、リージョンでデータセットをマスクする必要があります。3次元xarrayデータセットからeu_countries_gdfデータフレーム(NUTSシェープファイル)、経度、緯度を渡し、eu_countries_gdfデータフレームのインデックス列に従って、国ごとに一意の整数値を含む一意の番号で特定の地域をマークします。

次に、マスクされた領域を選択し、特定の月(この場合は6月)のみを選択します。

import regionmask
# Create mask of multiple regions from shapefile
eu_mask = regionmask.mask_3D_geopandas(
        eu_countries_gdf,
        multi_ds.longitude,
        multi_ds.latitude,
        drop=True,
        numbers="index"
    )
# Apply mask on our dataset
multi_ds = multi_ds.where(eu_mask)
# Get indexes of months since we have January and June
# Select January (1 : January, 6 : June etc...)
month_idxs = multi_ds.groupby('valid_time.month').groups
ds_1month = multi_ds.isel(valid_time=month_idxs.get(6))
ds_1month

      
                

grouped_ds = multi_ds.groupby("valid_time.year").mean(["latitude","longitude"])
grouped_ds

次に、結果をpandasデータフレームに変換し、いくつかの変換を実行して、結果を視覚化します。

  1. dfeu_countries_gdfと結合して、国の文字列IDであるNUTS_IDに関する情報を遡及的に取得できるようにします
  2. このデータセットではvalid_timeと同じであるため、列の時間を削除します
  3. valid_timeをインデックスとして設定し、irに従ってデータを並べ替えます
  4. NUTS_IDのピボットデータフレーム
  5. t2mパラメータを選択し、軸1の名前を「国」に変更します
  6. 日時 インデックスをの整数に変換します
  7. インデックスの名前を「Year」に変更
  8. df = grouped_ds.to_dataframe()
    df = df.reset_index("valid_time").join(eu_countries_gdf.set_index("index")[["NUTS_ID"]])
    df = df.drop(columns=["time"])
    df = df.set_index("valid_time").sort_index()
    df = df.pivot(columns=["NUTS_ID"])
    df = df["t2m"] # Select temperature (this is useful if we have more parameters in grid, eg. wind, humidity)
    df = df.rename_axis("Country", axis=1)
    df.index = [value.year for value in df.index]
    df.index.name = 'Year'
    df.head()
    
           
                    

    import geopandas as gpd
    import xarray as xr
    import regionmask
    import glob
    import pandas as pd
    def prepare_data(multi_ds, month : int = 1) -> pd.DataFrame:
        
        month_idxs = multi_ds.groupby('valid_time.month').groups
        multi_ds = multi_ds.isel(valid_time=month_idxs.get(month))
        multi_ds = multi_ds.groupby("valid_time.year").mean(["latitude","longitude"])
        df = multi_ds.to_dataframe()
        df = df.reset_index("valid_time").join(eu_countries_gdf.set_index("index")[["NUTS_ID"]])
        df["valid_time"] = df["valid_time"].apply(lambda x: x.year)
        df = df.drop(columns=["time"])
        df = df.set_index("valid_time").sort_index()
        df = df.pivot(columns=["NUTS_ID"])
        df = df["t2m"]
        df = df.rename_axis("Country", axis=1)
        df.index.name = 'Year'
        return df
    # Get list of files with anomaly data using glob
    grib_files = glob.glob("data/ecv/temperature/*anomaly*.grib")
    # Open multifile dataset with xarray, concat them on "time" dimension
    with xr.open_mfdataset(grib_files, combine="nested",concat_dim = "valid_time", parallel=True) as data:
        multi_ds = data
        
    # Open shapefile with geopandas
    eu_countries_shp = "data/shapefiles/NUTS_RG_01M_2021_4326.shp"
    eu_countries_gdf = gpd.read_file(eu_countries_shp)
    eu_countries_gdf = eu_countries_gdf[eu_countries_gdf["LEVL_CODE"] == 0]
    eu_countries_gdf = eu_countries_gdf.reset_index()
    # Create mask of multiple regions from shapefile
    eu_mask = regionmask.mask_3D_geopandas(
            eu_countries_gdf,
            multi_ds.longitude,
            multi_ds.latitude,
            drop=True,
            numbers="index"
        )
    # Apply mask on our dataset
    multi_ds = multi_ds.where(eu_mask)
    january_df = prepare_data(multi_ds, month=1)
    june_df = prepare_data(multi_ds, month=6)
    january_df.head()
    
           
                    

import matplotlib.pyplot as plt
import seaborn as sns
grid_kws = {"width_ratios": (.95, .03), "hspace": .1}
f, (ax, cbar_ax) = plt.subplots(1,2,figsize=(16,12), gridspec_kw=grid_kws)
ax = sns.heatmap(june_df.T, ax=ax,
                 cbar_ax=cbar_ax,
                 square=True, cmap='RdBu_r', vmin=-5, vmax=5, center=0, annot=True, fmt=".0f",
                 cbar_kws={"orientation": "vertical"})
ax.set_title("Monthly Mean Temperature Anomaly in June in Selected European Countries (°C)", size="18", pad=20)

ほとんどのヨーロッパ諸国で気温上昇の傾向をはっきりと見ることができます(著者による画像)

結果について考えてみましょう。横軸は昇順の年数を表しています。縦軸は、NUTS_IDでラベル付けされた国を表します。私たちのカラースケールは温度異常を表しており、負の値は青みがかった色で表示され、正の値は赤みがかった色で表示されます。

上の画像から、ヨーロッパのほとんどの国で、6月は10年から10年にかけて暖かくなっていると断言できます。1980年から1995年までの数年間、ほとんどの国の気温は通常より低くなりましたが、ここ2、3年は反対側で、気候平均よりもはるかに暖かくなりました(この場合、気候平均は、グリッド内の各ポイントの平均気温です。 1991年から2020年)。

1月の結果を見て、この冬の月の気温がどのように変化したかを判断しましょう。

非常に寒い1月の頻度が低下しています(著者による画像)

上の写真でわかるように、6月の気温ほど傾向ははっきりしていません。しかし、80年代後半には、2000年代よりも、1月の非常に寒い気温の頻度が高かったという結論に達することができます。

自分でデータセットを試してみて、調査結果を知らせてください。ボーナスとして、私はフランスの地域で6月に変化する気温も分析しました。データの解像度は0.25x0.25°であることに注意してください。これは、非常に小さな領域の分析には適していないことを意味します。

地域規模でも、フランスでは気温の上昇傾向が見られます。2003年の平均気温をはるかに上回っていることに注意してください(著者による画像)

この記事が気に入ったら、遠慮なく私のアカウントをフォローしてください。コメントは大歓迎です。ご質問にお答えさせていただきます。ご清聴ありがとうございました。今後ともよろしくお願いいたします。

接続を維持

  • このようなより多くの物語のためにミディアムで私に従ってください
  • LinkedInで接続する
  • CDSで気候変動に関するその他のデータセットを検索する

© Copyright 2021 - 2022 | hachiwiki.com | All Rights Reserved