Thực hiện thống kê vùng về dữ liệu khí hậu với Python

May 08 2022
Khí hậu của đất nước bạn thay đổi như thế nào so với các nước khác? Biến đổi khí hậu đang dần bắt đầu hình thành tương lai của tất cả chúng ta. Trong những thập kỷ tới, chúng ta có thể mong đợi nhiệt độ ấm hơn, hạn hán kéo dài, bão mạnh hơn và nhiều hơn nữa.

Khí hậu của đất nước bạn thay đổi như thế nào so với các nước khác?

Nhiệt độ tháng 6 đang tăng trên hầu hết các quốc gia ở Châu Âu (Ảnh của tác giả)

Biến đổi khí hậu đang dần bắt đầu hình thành tương lai của tất cả chúng ta. Trong những thập kỷ tới, chúng ta có thể mong đợi nhiệt độ ấm hơn, hạn hán kéo dài, bão mạnh hơn và nhiều hơn nữa.

Chúng tôi sẽ không thảo luận về biến đổi khí hậu trong bài viết này vì bạn có thể tìm thấy rất nhiều cuộc thảo luận về chủ đề này trên internet. Tuy nhiên, chúng tôi sẽ khám phá dữ liệu miễn phí từ Ủy ban Châu Âu và Copernicus có thể giúp chúng tôi hiểu nhiệt độ thay đổi như thế nào trong những năm qua.

Trong bài viết này, chúng tôi sẽ thực hiện thống kê vùng trên nhiều tập dữ liệu. Chúng tôi sẽ hình dung sự bất thường của nhiệt độ trung bình hàng tháng và tạo chuỗi thời gian trên nhiều quốc gia. Cuối cùng, chúng ta sẽ hình dung những phát hiện của mình bằng cách sử dụng biểu đồ bản đồ nhiệt.

TLDR: Nếu bạn không quan tâm đến cách lấy dữ liệu và môi trường thiết lập, hãy chuyển đến phần Thực hiện thống kê vùng và tận hưởng phần viết mã của bài viết.

Thiết lập môi trường

Thiết lập môi trường Python

Chúng ta sẽ bắt đầu với việc tạo môi trường python mới bằng Anaconda. Chúng tôi sẽ sử dụng nhiều thư viện không gian địa lý như xarray, regionmask hoặc geopandas. Chúng tôi cũng sẽ sử dụng cdsapi để kết nối với Kho dữ liệu khí hậu và lấy dữ liệu.

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

conda activate climate_change

Điều tiếp theo sẽ là thiết lập chính xác API CDS Python, mà chúng tôi sẽ sử dụng để tải xuống dữ liệu từ Climate Data Store. Hãy tiếp tục theo bài viết cuối cùng (và đầu tiên) của tôi .

Vui lòng tạo thư mục mới nơi bạn đặt tất cả dữ liệu và sổ ghi chép ipython, trong đó chúng tôi sẽ thực hiện tất cả các bước bên dưới. Tôi đã tạo sổ tay Clio_change.ipynb mà tôi sẽ chia sẻ ở cuối bài viết này.

Nhận dữ liệu mà chúng tôi sử dụng

Chúng tôi sẽ sử dụng lại kiến ​​thức từ phần trước của tôi và chúng tôi sẽ khám phá Kho dữ liệu khí hậu một lần nữa.

Nhận dữ liệu khí hậu

Hôm nay, chúng ta sẽ làm việc với Các biến khí hậu thiết yếu để đánh giá sự biến đổi khí hậu từ Copernicus. Tập dữ liệu bao gồm nhiều biến khí hậu như Nhiệt độ không khí bề mặt, Độ ẩm tương đối của không khí bề mặt hoặc Lượng mưa. Sau đó, chúng tôi có lựa chọn Bất thường (Khí hậu trừ đi giá trị trung bình hàng tháng), Khí hậu hoặc loại sản phẩm trung bình hàng tháng .

Do lượng dữ liệu lớn, chúng tôi sẽ so sánh nhiệt độ của tháng 1 và tháng 6 trong suốt 44 năm từ năm 1979 đến năm 2022 (lưu ý rằng chúng tôi sẽ chỉ có 43 bộ dữ liệu cho tháng 6, vì bài báo này được viết vào tháng 5 năm 2022). Cuối cùng thì chúng ta sẽ có 87 tệp.

Yêu cầu API python của chúng tôi sẽ giống như sau:

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')

Sau khi dữ liệu được tải xuống, bạn nên giải nén nó vào thư mục bạn chọn. Do đó, bạn sẽ có 87 tệp grib được đặt trong thư mục mong muốn của mình.

Tạo / Lấy các tệp hình dạng

Để tính toán số liệu thống kê theo vùng, chúng tôi muốn tính toán số liệu thống kê trên một khu vực cụ thể. Trong trường hợp của chúng tôi, các khu vực cụ thể là các nước Châu Âu. Chúng tôi sẽ sử dụng NUTS shapefile - một tệp chứa thông tin không gian địa lý về các quốc gia thuộc Liên minh Châu Âu.

Chúng tôi sẽ sử dụng shapefile, tệp đã được tạo bởi Ủy ban Châu Âu. Tuy nhiên, nếu bạn đã quen với hệ thống GIS, bạn có thể tạo một hệ thống của riêng mình.

Bộ dữ liệu shapefile với mô tả của chúng tôi có sẵn để tải xuống tại đây . Chúng tôi sẽ sử dụng cấu hình này và tải xuống.

Cấu hình để tải xuống shapefile được sử dụng trong bài viết này

Lưu nó vào vị trí ưa thích, tôi thực sự khuyên nó nên là cùng một thư mục với thư mục mà chúng tôi lưu trữ sổ ghi chép python.

NUTS là một hệ thống phân cấp được chia thành 3 cấp độ. NHÓM 1: các vùng kinh tế - xã hội chính, NÓI 2: các vùng cơ bản để áp dụng các chính sách vùng, NHÓM 3: các vùng nhỏ để chẩn đoán cụ thể.

Chúng tôi sẽ sử dụng NUTS1 và trong ví dụ thưởng, chúng tôi sẽ làm việc với NUTS2.

Đây là cách shapefile NUTS của chúng tôi trông như thế nào trong QGIS, thêm một số ký hiệu cơ bản để phân biệt giữa các quốc gia:

Bộ dữ liệu NUTS được hiển thị trong QGIS (Hình ảnh của Tác giả)

Thực hiện thống kê vùng

Vì chúng tôi đã thu thập tất cả dữ liệu mà chúng tôi cần, chúng tôi có thể tiếp tục với bước thứ hai đến bước cuối cùng - Thực hiện Thống kê Vùng. Những gì chúng tôi muốn đạt được là tính toán nhiệt độ trung bình cho mọi phần tử (quốc gia) trong tệp hình dạng của chúng tôi trong suốt 44 năm cho tháng Giêng và 43 năm cho tháng Sáu, tương ứng. Chúng tôi sẽ chỉ thực hiện tính toán cho loại sản phẩm Anomaly (bạn có thể tự tính toán Giá trị trung bình hàng tháng sau khi đọc bài viết này). Chúng tôi sẽ trích xuất thông tin này vào Pandas DataFrame, mà chúng tôi sẽ sử dụng để hiển thị kết quả cuối cùng.

Như một phần thưởng, chúng tôi sẽ chọn một quốc gia và thực hiện thống kê theo vùng trên các vùng khác nhau trong quốc gia.

Hãy đi sâu vào phần viết mã của bài viết này.

Khám phá tập dữ liệu khí hậu

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

      
                

Tuy nhiên, điều này có thể được thực hiện dễ dàng với hàm xarray open_mfdataset () , nơi chúng ta có thể chỉ định danh sách các tệp để mở và cũng chỉ định kích thước lồng nhau và kết hợp. Trong trường hợp của chúng tôi, chúng tôi muốn nối dữ liệu qua valid_time , mà chúng tôi có thể thấy dưới tọa độ trong tập dữ liệu ở trên. Đặt các đối số phù hợp cho open_mfdataset () sẽ giúp chúng tôi mở tất cả các tệp mà chúng tôi cần, thẳng đến tập dữ liệu 3 chiều. Ngoài ra, tất cả các thao tác mở có thể được thực hiện song song trên nhiều bộ xử lý.

# 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]

Để tính giá trị trung bình theo vùng qua các vùng, trước tiên chúng ta phải che tập dữ liệu với các vùng của chúng ta, sử dụng mặt nạ vùng thư viện. Chúng tôi sẽ chuyển khung dữ liệu eu_countries_gdf (NUTS shapefile), kinh độ và vĩ độ từ tập dữ liệu xarray 3 chiều của chúng tôi và chúng tôi sẽ đánh dấu các vùng cụ thể của mình bằng số duy nhất theo cột chỉ mục trong khung dữ liệu eu_countries_gdf , chứa các giá trị số nguyên duy nhất cho mọi quốc gia.

Sau đó, chúng tôi sẽ chọn các khu vực được che và chỉ chọn tháng cụ thể (trong trường hợp này là tháng 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

Bây giờ, hãy chuyển đổi kết quả của chúng tôi sang khung dữ liệu gấu trúc, thực hiện một số phép biến đổi và trực quan hóa những phát hiện của chúng tôi:

  1. Tham gia df với eu_countries_gdf để chúng tôi có thể truy xuất thông tin về NUTS_ID, là chuỗi ID của các quốc gia
  2. Thả thời gian của cột , vì trong tập dữ liệu này, nó giống như valid_time
  3. Đặt valid_time làm chỉ mục và sắp xếp dữ liệu theo ir
  4. Tổng hợp khung dữ liệu trên cột NUTS_ID
  5. Chọn tham số t2m , đổi tên trục 1 thành “Quốc gia”
  6. Chuyển đổi chỉ mục ngày giờ thành số nguyên của năm
  7. Đổi tên chỉ mục thành “Năm”
  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)

Chúng ta có thể thấy rõ xu hướng tăng nhiệt độ ở hầu hết các nước Châu Âu (Ảnh của tác giả)

Hãy suy nghĩ về kết quả. Trục hoành của chúng tôi đại diện cho các năm theo hướng tăng dần. Trục tung của chúng tôi đại diện cho các quốc gia được gắn nhãn NUTS_ID. Thang màu của chúng tôi thể hiện sự bất thường về nhiệt độ với các giá trị âm có màu hơi xanh, trong khi các giá trị dương được hiển thị bằng màu hơi đỏ.

Từ hình ảnh trên, chúng ta chắc chắn có thể nói rằng ở hầu hết các quốc gia Châu Âu, nước trái cây đang trở nên ấm hơn từ thập kỷ này sang thập kỷ khác. Trong khi trong những năm từ 1980 đến 1995, hầu hết các quốc gia đều có nhiệt độ dưới mức bình thường, thì vài năm qua lại ở phía ngược lại, ấm hơn nhiều so với mức trung bình của khí hậu (trong trường hợp của chúng tôi, trung bình khí hậu là nhiệt độ trung bình cho mỗi điểm trong lưới được tính trên năm 1991–2020).

Hãy cùng xem kết quả của tháng 1 để xác định, nhiệt độ đã thay đổi như thế nào trong tháng mùa đông này.

Tần suất các tháng Giêng rất lạnh đang giảm (Ảnh của tác giả)

Như chúng ta có thể thấy trên hình trên, xu hướng không rõ ràng như vậy đối với nhiệt độ vào tháng Sáu. Tuy nhiên, chúng ta có thể kết luận rằng vào cuối những năm 80, tần suất nhiệt độ rất lạnh trong tháng 1 cao hơn so với những năm 2000.

Hãy thử tự xử lý tập dữ liệu và cho tôi biết những phát hiện của bạn. Phần thưởng, tôi cũng đã phân tích sự thay đổi nhiệt độ trong tháng 6 ở các vùng của Pháp. Lưu ý rằng độ phân giải của dữ liệu là 0,25x0,25 °, có nghĩa là nó không phù hợp để phân tích các vùng rất nhỏ.

Ngay cả trên quy mô khu vực, chúng ta có thể thấy xu hướng nhiệt độ tăng trên Pháp. Lưu ý nhiệt độ cao hơn nhiều so với trung bình trong năm 2003 (Ảnh của tác giả)

Nếu bạn thích bài viết này, đừng ngần ngại theo dõi tài khoản của tôi. Mọi nhận xét đều được hoan nghênh và tôi rất vui được trả lời bất kỳ câu hỏi nào của bạn. Cảm ơn bạn đã chú ý và tôi mong muốn khám phá thêm dữ liệu với bạn.

Giữ liên lạc

  • Theo dõi tôi trên Medium để biết thêm nhiều câu chuyện như thế này
  • Kết nối trên LinkedIn
  • Tìm thêm bộ dữ liệu về Biến đổi khí hậu trên CDS

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