[파이썬] 전구 자료 위경도 평균해보자! 추세선까지 그려보자!

안녕하세요? 째까입니다.

 

이번 파이썬 과정은 전구 자료를 가지고 x축은 시간, y축은 기온 자료로 시계열 그래프를 만들어 보는 것입니다!

 

자료는 NCEP-NCAR Reanalysis 1 자료를 사용했구요? monthly 자료입니다!

 

밑에서 다운 받으실 수 있슴다.

https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html

: NOAA Physical Sciences Laboratory

 

 

저는 이걸 선택했습니다~

열어보시면

이렇게 되어있는데요? level은 선택하시고.. lat, lon은 평균화합시다! 이번 글의 핵심이죠

 

먼저 자료를 읽어볼까요?

필요한 라이브러리를 import 합니다.

from warnings import filterwarnings
filterwarnings(action='ignore', category=DeprecationWarning, message='`np.bool` is a deprecated alias')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams.update({"font.size": 18})
import matplotlib.patches as patches

import matplotlib.font_manager as fm
plt.rcParams['font.family'] = "Times New Roman"
plt.rcParams["axes.unicode_minus"] = False

import os
import xarray as xr

from scipy import signal

그 다음 파일을 열죠

def open_file(temp_name, sst_name) :
    tmp_daio = os.path.join(os.getcwd(), "DAIO")
    tmp_dir = os.listdir(tmp_daio)

    temp_file = [file for file in tmp_dir if file.startswith(temp_name) and "_6190_" in file and "pres" in file]
    temp_nc_path = os.path.join(tmp_daio, *temp_file)
    temp_nc = xr.open_dataset(temp_nc_path)
    
   
    return temp_nc

파일은 temp_file로 조건에 맞는걸 필터한 후 temp_nc에 할당한 후에! xarray로 읽어서 return했슴당

 

이제 lat, lon을 평균할껀데 그냥 mean하면 안된다네요?

lat의 사이즈가 전구 모두 같지 않기 때문에!! 아래 방법을 써야한다고 합니다 ㅎㅎ.. 으휴

def weights_var(ds, syy, eyy, season, lwhpa = None, hihpa = None, sst = "False"):
    if sst == "False" :
        ss_ds = select_yy_mm(ds, syy, eyy, season)
        ss_ds = ss_ds.groupby("time.year").mean("time")
        ss_ds = ss_ds.sel(level = hihpa) - ss_ds.sel(level = lwhpa)

    elif sst != "False" :
        ss_ds = select_yy_mm(ds, syy, eyy, season)
        ss_ds = ss_ds.groupby("time.year").mean("time")


    weights = np.cos(np.deg2rad(ss_ds.lat))
    weights.name = "weights"
    tmp_weighted = ss_ds.weighted(weights)
    ss_mean = tmp_weighted.mean(("lon", "lat"))

    return ss_mean

저는 SST도 사용해서 sst 조건문이 들어가있습니다ㅎㅎ... 무시해주시고

select_yy_mm는 제가 원하는 기간과 월만 선택한거니 무시해주십쇼!

저는 무시하지마십쇼!

 

무튼 바로 위 코드에서 중요한건 weights부분인데요. nc파일 중에 lat을 선택한 후 weight하는 것입니다!

lat weight 후 lon과 lat을 평균합니당 그럼 이제 time만 남겠져?(level은 이미 선택했습니다!)

 

 

 

그 다음에 추세선을 계산합니당

def cal_time_poly(nc, diff = "False") :
    x = nc["year"]
    x2 = np.arange(1, len(x)+1)
    x3 = np.arange(len(x))
    y = nc["air"]

    z = np.polyfit(x2, y,1).round(2)
    pred = np.poly1d(z)

    if diff == "False" :
        return x2, y, pred

이런식으로 x는 시간 y는 어떤 level의 기온입니다!

회귀분석은 numpy의 polyfit을 이용합니다

polyfit은(x, y, 차수)가 기본이고요? 저는 1차 방정식이용할꺼니까 1을 넣었슴다

그러면 z에 y = ax+b에서 a랑 b가 return되고, a, b를 따로 따로 받고싶으시면 

a, b = np.polyfit(x2, y, 1)이런식으로 쓰세욥!

 

그 다음 poly1d에 z를 넣어서 y = ax + b를 가지고있고(pred)

return합니다! 나중에 pred에 값(개수)을 넣어서 추세 그림을 그립니당

---

자 이제 위에 쓴 함수들을 call하고 그림 그립시다!

다왔어요 여러분

 

def check_poly_plot(temp_nc, sst_nc):
    global lwhpa, hihpa, wnt, smr, daou_path
    wnt_temp_ds = weights_var(temp_nc, syy, eyy, wnt, lwhpa, hihpa)
    wnt_time = wnt_temp_ds["year"]
    wnt_temp = wnt_temp_ds["air"]
    

    wnt_x2, wnt_orig_temp, wnt_orig_reg = cal_time_poly(wnt_temp_ds)
    

    fig1, ax1 = plt.subplots(figsize = (14, 8))
    ax1.plot(wnt_time, wnt_orig_temp, color = "r", label = "Origin")
    ax1.plot(wnt_time, wnt_orig_reg(wnt_x2), color = "r", linestyle = "--")

    ax1.set_xlabel("Year")
    ax1.set_ylabel("Temperature")
    ax1.legend()

    plt.tight_layout()
    plt.show()

 

이제 보이시죠?

weights_var로 위경도 평균하고?

cal_time_poly로 추세 분석하고?

 

그림 그리는 겁니다 간단하죠?

 

그림 그리면?

짜잔!

완성!

 

다음에 또 찾아올께욥