개화모델 구현 및 시각화 프로젝트

개화모델 구현 및 시각화 프로젝트

프로젝트 개요

이 프로젝트에서는 사과와 배의 개화 시점을 예측하기 위한 여러 개화 모델을 파이썬 코드로 구현하고, 이를 GUI 환경에서 시각화하였다. 식물의 개화는 농업에서 방제 계획 수립 등에 있어 중요한 정보이다. 개화 모델은 온도, 일장 등 환경 요인과 식물의 생리적 반응을 수학적으로 연결하여 개화 시기를 예측하는 데 활용한다. 이 모델들을 사용하여 개화 시점을 예측하고 이를 시각적으로 확인할 수 있도록 하였다.

개화 모델 선정

1. DVR(Degree-Day Vegetation Growth Model)

DVR 모델은 특정 온도 이상에서 식물이 성장하는 데 필요한 ‘degree-day’를 계산하여 개화 시점을 예측한다. 누적된 발달 속도가 특정 임계값에 도달하면 개화가 일어난다고 가정한다.

• 입력자료: 일평균기온
• 발육속도(development rate, DVR): 과수가 일평균 기온에 영향을 받아 만개기에 다가가는 발육속도
• 발육단계(development stage, DVS): 발육속도에 의해 도달되는 발육단계
• 예상 만개기: 누적되는 DVR 값이 100에 도달하면 해당 날짜를 예상 만개기로 출력

# -------DVR 모델
def dvr_model(df):
    A = 107.94
    B = 0.9

    # 발육속도 > 5°C 이상의 온도만 계산
    df['DVRi'] = df['temp'].apply(lambda x: 100 / (A * (B ** x)) if x >= 5 else 0)

    # 발육단계 > 누적발육속도(∑DVRi. DVS)
    df['DVS'] = df.groupby('year')['DVRi'].cumsum()

    # 예상만개기 > 발육단계 == 100
    expected_full_bloom = df[df['DVS'] >= 100].groupby('year').first()

    dvr_date = expected_full_bloom['date'].reset_index()
    dvr_date = dvr_date.rename(columns={"date" : "DVR"})
    return dvr_date

2. mDVR(Modified DVR)

DVR 모델을 개선한 버전으로, 더 정확한 예측을 위해 추가 요인들을 고려한다. 예를 들어, 춘화 요구도 등 품종 특성을 반영하여 모델을 보완한다.

• 입력자료: Sugiura(1999)의 방법에 따라 전날의 최고기온과 당일 최저기온으로부터 시간별 기온을 추정

# -------정시의 기온 예측
def predict_time_temp(df):
    for hour in range(0, 24):
        if 0 <= hour <= 3:
            df[f'temp_{hour}시'] = (df['max_temp'].shift(1) - df['min_temp']) * (np.sin((4 - hour) * 3.14 / 30) ** 2) + df[
                'min_temp']

        elif 4 <= hour <= 13:
            df[f'temp_{hour}시'] = (df['max_temp'] - df['min_temp']) * (np.sin((hour - 4) * 3.14 / 18) ** 2) + df['min_temp']

        elif 14 <= hour <= 23:
            df[f'temp_{hour}시'] = (df['max_temp'] - df['min_temp'].shift(-1)) * (np.sin((28 - hour) * 3.14 / 30) ** 2) + df[
                'min_temp'].shift(-1)
    return df

• 내생휴면 발육속도(DVR1): 휴면타파에 필요한 0∼6℃ 범위의 적산 시간을 온도범위에 따른 발육속도(DVR1)으로 환산

# -------시간별 DVR1
def get_dvr1(x):
    dvr = None
    if x < -6:
        dvr = 0

    elif -6 <= x < 0:
        dvr = 1.333 * (10 ** (-3)) + 2.222 * (10 ** (-4)) * x

    elif 0 <= x < 6:
        dvr = 1.333 * (10 ** (-3))

    elif 6 <= x < 9:
        dvr = 2.276 * (10 ** (-3)) - 1.571 * (10 ** (-4)) * x

    elif 9 <= x < 12:
        dvr = 3.448 * (10 ** (-3)) - 2.784 * (10 ** (-4)) * x

    elif 12 <= x:
        dvr = 0

    return dvr

• 강제휴면 발육속도(DVR2): 저온감응기 이후 만개기에 도달할 때까지 필요한 발육속도(DVR2)

  # -------시간별 DVR2
  def get_dvr2(x):
    dvr2 = None
  
    if x < 20:
        a = 35.27 - 12094 * ((x + 273) ** (-1))
        dvr2 =  np.exp(a)
    elif x >= 20:
        a = 5.82 - 3474 * ((x + 273) ** (-1))
        dvr2 =  np.exp(a)
  
    elif x < 0:
        dvr2 = 0
  
    return dvr2
  

• 내생휴면 해제: 누적발육지수(∑DVR1)이 1에 도달하면 내생휴면해제
• 저온감응기: 내생휴면과 강제휴면이 겹치는 시기로, 누적발육지수(∑DVR1)이 2에 도달하면 종료
• 예상만개기: 누적되는 강제휴면 발육속도(∑DVR2) 값이 0.9593에 도달하면 해당날짜를 예상만개기로 출력

def modified_dvr_model(df):
    df['season_year'] = df['year']
    df.loc[(df["month"] >= 10), "season_year"] = df["year"] + 1
    df = df[~df['month'].isin([7, 8, 9])]
    df = df.sort_values(by=['year', 'date'], ascending=[True, True])
    df = predict_time_dvr(df)

    all_date = pd.DataFrame()

    # 내생휴면해제종료일 2월 15일로 설정
    for year in df['season_year'].unique():
        year_df = df[df['season_year'] == year]

        # col_dvr1 = [col for col in year_df.columns if 'DVR1' in col]
        # year_df['DVR1'] = year_df[col_dvr1].sum(axis=1)
        # year_df['DVS1'] = year_df['DVR1'].cumsum()
        # year_df = year_df[year_df['DVS1'] >= 1] # 내생휴면 해제

        # year_df = year_df[year_df['DVS1'] >= 2] # 저온감응기(내생휴면과 강제휴면이 겹치는 시기) 종료 이후

        year_df = year_df[(year_df['date'] >= pd.to_datetime(f'{year}-02-15'))] # 내생휴면 해제 종료일 설정
        # col_dvr2 = [col for col in year_df.columns if 'DVR2' in col]
        # year_df['DVR2'] = year_df[col_dvr2].sum(axis=1) # 만개기에 도달할 때까지 필요한 발육속도
        # year_df['DVS2'] = year_df['DVR2'].cumsum() # 발육속도 누적
        year_df['DVS2'] = year_df['dvr2'].cumsum() # 발육속도 누적

        year_date = year_df[year_df['DVS2'] >= 0.9593].head(1) # 해당 연도의 예측 만개기 출력
        all_date = pd.concat([all_date, year_date])

    mdvr_date = all_date[['season_year', 'date']].reset_index(drop=True)
    mdvr_date = mdvr_date.rename(columns={"date" : "mDVR",
                                          "season_year": "year"})

    return mdvr_date

3. CD(Chill Days)

저온 요구도를 충족시키는 데 필요한 저온 일수를 계산한다. 일정 온도 이하의 날짜 수를 누적하여 개화에 필요한 저온 요구도를 산출한다.

• 입력자료: 매일의 최고기온과 최저기온을 이용해 기준온도로부터 유효한 온도범위에 따라 가중치를 달리해 내생휴면 해제 이전은 냉각량(chill unit, CU), 그 이후는 가온량(heat unit, HU)로 표현

def add_chill_heat(row):
    tc = 5.4
    max = row['max_temp']
    min = row['min_temp']
    avg = row['temp']

    chill = 0
    anti_chill = 0

    if 0 <= tc <= min <= max:
        chill = 0
        anti_chill = avg - tc
    elif 0 <= min <= tc <= max:
        chill = -((avg - min) - (max - tc) / 2)
        anti_chill = (max - tc) / 2
    elif 0 <= min <= max <= tc:
        chill = -(avg - min)
        anti_chill = 0
    elif min < 0 <= max <= tc:
        chill = (max / (max - min)) * (max / 2)
        anti_chill = 0
    elif min < 0 < tc < max:
        chill = -(((max / (max - min)) * (max / 2)) - ((max - tc) / 2))
        anti_chill = (max - tc) / 2


    row['chill_unit'] = chill
    row['heat_unit'] = anti_chill

    return row

• 내생휴면해제: 휴면에 진입 후 누적 냉각량이 정해진 저온요구도(chilling requirement, Cr)에 도달
• 강제휴면타파(발아): 내생휴면해제 이후 누적 가온량이 저온요구량에 도달
• 예상만개기: 신고 배의 고온요구도(heating requirement, Hr)만큼 추가적인 가온량 누적

def cd_model(df):
    cr = -86.4 # 저온요구도
    hr = 272 # 고온요구도

    df['date'] = pd.to_datetime(df['date'])
    df['season_year'] = df['year']
    df.loc[(df["month"] >= 10), "season_year"] = df["year"] + 1
    df = df[~df['month'].isin([7, 8, 9])]
    df = df.sort_values(by=['year', 'date'], ascending=[True, True])

    # chill days 계산
    df = df.apply(add_chill_heat, axis=1)

    all_date = pd.DataFrame()
    for year in df['season_year'].unique():
        year_df = df[df['season_year'] == year]

        # 내생 휴면 해제
        # year_df['cumsum_chill_unit'] = year_df['chill_unit'].cumsum()
        # year_df = year_df[year_df['cumsum_chill_unit'] >= cr]
        year_df = year_df[(year_df['date'] >= pd.to_datetime(f'{year}-02-15'))]

        year_df['cumsum_heat_unit'] = year_df['heat_unit'].cumsum()
        year_df = year_df[year_df['cumsum_heat_unit'] >= cr]

        # 만개예상일
        # year_df['cumsum_heat_unit'] = year_df['heat_unit'].cumsum()
        year_date = year_df[year_df['cumsum_heat_unit'] >= hr].head(1)

        all_date = pd.concat([all_date, year_date])

    cd_date = all_date[['year', 'date']].reset_index(drop=True)
    cd_date = cd_date.rename(columns={"date": "CD"})

    return cd_date

4. GDH(Growth Degree Hours)

식물 생장에 필요한 열량을 시간 단위로 누적하여 계산한다. 기준 온도 이상의 온도에서 1시간당 1GDH가 누적된다. 누적된 GDH가 특정 임계값에 도달하면 개화가 일어난다고 예측한다.

• 입력자료: Sugiura(1999)의 방법에 따라 전날의 최고기온과 당일 최저기온으로부터 시간별 기온을 추정(Table2)
• 휴면타파시점 이후 시간별 기온을 GDH로 변환∙대입한 뒤 개화추정모델의 기준치(25, 30, 50%; 각 GDH 1395, 1674, 2790)에 도달한 시점으로 예상되는 개화시기를 추정
• 후지 사과 개화추정모델: GDH5579

def cal_gdh_value(x):
    gdh = 0
    if 4 < x < 25:
        gdh = 10.5 * (1 + np.cos(np.pi + np.pi * ((x - 4) / 21)))
    elif 25 < x < 36:
        gdh = 21 * (1 + np.cos((np.pi / 2) + (np.pi / 2) * ((x - 4) / 25)))
    return gdh

def predict_time_gdh(df):
    id_vars = [col for col in df.columns if not 'temp_' in col]

    for hour in range(0, 24):
        df[f'gdh_{hour}시'] = df[f'temp_{hour}시'].apply(cal_gdh_value)

    col_gdh = [col for col in df.columns if 'gdh_' in col]
    df = df.melt(id_vars=id_vars, value_vars=col_gdh, var_name="time", value_name='gdh')
    df['time'] = df['time'].apply(lambda x: re.findall(r'\d+', x)[0]).astype(int)
    df = df.sort_values(by=['date', 'time'], ascending=[True, True]).reset_index().drop(columns='index')

    return df

def gdh_model(df):
    # after data 필요
    gdh_values = [1395, 1674, 2790, 5579]

    after_df = pd.DataFrame()

    for year in df['year'].unique():
        # year = row['year']
        year_df = df[df['year'] == year]
        year_df = year_df[year_df['date'] >= pd.to_datetime(f'{year}-02-20')]

        after_df = pd.concat([after_df, year_df], ignore_index=True)

    # after = pd.concat(after_df)
    after = predict_time_gdh(after_df)
    after['cumsum_gdh'] = after.groupby('year')['gdh'].cumsum()

    dates_list = []

    for gdh_value in gdh_values:
        dormancy_dates = after[after['cumsum_gdh'] >= gdh_value].groupby('year')['date'].first().reset_index()
        dormancy_dates = dormancy_dates.rename(columns={'date': f'gdh_{gdh_value}'})
        dates_list.append(dormancy_dates)

    gdh_date = pd.concat(dates_list, axis=1)
    gdh_date = gdh_date.loc[:, ~gdh_date.columns.duplicated()]
    return gdh_date

5. Utah

온도에 따른 저온 요구도 충족 효과를 가중치로 부여하여 계산한다. 각 온도 구간별로 다른 가중치(Chill Unit)을 적용한다. 누적된 Chill Unit이 품종별 요구도에 도달하면 휴면이 타파되고 개화로 진행된다고 본다.

• 온대과수의 휴면: 식물체 내부적으로 생리적 변화가 있으나 외관상 생장이 멈춘 상태
• 휴면의 구분: 외재휴면, 내재휴면, 환경휴면 단계
• 내재휴면: 계절적으로 식물 생장에 불리한 조건을 극복하기 위한 수체의 내부적인 생리반응으로 휴면타파를 위해서는 저온 축적이 필요
• 저온요구도: 과종, 품종에 따라 다르며 5~7.2℃ 내외에서 효과
• 겨울철 불충분한 저온의 축적으로 발아와 개화가 지연되거나 불량해지면 과실 생산량에 영향
• 사과: 500~800CU, 배: 1,000~1,200CU
• 저온요구도는 Utah모델에서 제시한 유효 온도별 가중치를 부여해 Chill Unit(CU)을 계산

def add_chill_unit(x):
    if x <= 1.4:
        return 0
    elif 1.5 <= x <= 2.4:
        return 0.5
    elif 2.5 <= x <= 9.1:
        return 1.0
    elif 9.2 <= x <= 12.4:
        return 0.5
    elif 12.5 <= x <= 15.9:
        return 0
    elif 16.0 <= x <= 18.0:
        return -0.5
    elif 18 <= x:
        return -1.0

def predict_chill_unit(df):
    for hour in range(0, 24):
        df[f'chill_{hour}시'] = df[f'temp_{hour}시'].apply(add_chill_unit)

    col_chill = [col for col in df.columns if 'chill_' in col]
    df['date_cumsum_chill'] = df[col_chill].sum(axis=1) # 일별 저온요구도
    df = df.drop(columns=col_chill)
    return df

• 휴면개시일: 저온요구도가 처음으로 0이상이 된 날
• 휴면타파일: 휴면개시일 이후 저온축적량이 저온요구도에 도달하면 휴면타파일로 출력

def utah_model(df,cu_MAX):
    df = predict_chill_unit(df)

    # 휴면개시일
    df['season_year'] = df["year"]
    df.loc[df["month"] > 9, "season_year"] = df["year"] + 1
    df = df.sort_values(by=['year', 'season_year', 'month'], ascending=[True, False, True])
    dormancy_dates = df[df['date_cumsum_chill'] >= 0].groupby('season_year')['date'].first().reset_index()

    start_dormancy = []
    for idx, row in dormancy_dates.iterrows():
        season_year = row['season_year']
        date = row['date']
        year_df = df[df['season_year'] == season_year]
        year_df = year_df[(year_df['date'] >= pd.to_datetime(date))]
        start_dormancy.append(year_df)

    after = pd.concat(start_dormancy)
    after = predict_chill_unit(after)
    after['chill_unit_cumsum'] = after.groupby('season_year')['date_cumsum_chill'].cumsum()

    stop_dormancy_dates = after[after['chill_unit_cumsum'] >= cu_MAX].groupby('season_year')['date'].first().reset_index()
    dormancy_dates = pd.merge(dormancy_dates, stop_dormancy_dates, on='season_year', how='inner')
    dormancy_dates = dormancy_dates.rename(columns={'season_year': 'year',
                                                    'date_x': 'start_dormancy_date',
                                                    'date_y': 'stop_dormancy_date', })

    return dormancy_dates

과수별 개화 모델 적용 및 내생휴면해제 종료일 선정

내생휴면 해제 후 냉각량 적산에서 가온량 적산으로 전환하는 순차휴면모델은 따뜻한 지역에서 겨울철 충분하지 않은 저온으로 인해 휴면이 불완전한 상태로 유지가 되고, 이 후 현실성이 떨어지는 개화 예측일을 도출할 수 있기 때문에 휴면 과정을 포함하는 mDVR 모델과 CD 모델, GDH은 내생휴면해제 종료일을 추가로 설정하였다.

자료 수집

1. 기상청(종관기상) 관측 지점 자료 수집

기상청에서는 종관기상관측 장비로 관측한 일 기상자료를 조회하는 서비스를 제공하고 있다. 관측 지점 번호와, 시작일, 종료일 값을 지정하여 요청하면 평균 상대습도, 평균 온도, 총 강우량 등의 데이터를 얻을 수 있다.

import json
import requests
import pandas as pd
from urllib.parse import quote_plus, urlencode
url = 'http://apis.data.go.kr/1360000/AsosDalyInfoService/getWthrDataList'
servicekey = 'API_KEY'

headers = {'User-Agent': 'Mozilla/5.0 (Windows NT 6.3; Win64; x64)'
                     'AppleWebKit/537.36 (KHTML, like Gecko) Chrome/63.0.3239.132'
                     'Safari/537.36'}

params = f'?{quote_plus("ServiceKey")}={servicekey}&' + urlencode({
quote_plus("pageNo"): "1", 
quote_plus("numOfRows"): "720",
quote_plus("dataType"): "JSON",
quote_plus("dataCd"): "ASOS",
quote_plus("dateCd"): "DAY",
quote_plus("startDt"): start_date,
quote_plus("endDt"): end_date,
quote_plus("stnIds"): station_id
})

result = requests.get(url + params, headers=headers)

js = json.loads(result.content)
df = pd.DataFrame(js['response']['body']['items']['item'])

GUI 시각화

일평균 기온 변화 추이: 평년 vs. 선택연도

저온 축적량과 휴면타파시기

만개일 예측 결과

만개월의 평균 기온과 만개일 변화 추이

만개일의 일평균 기온과 만개일 변화 추이

전국 만개일과 만개월 평균 기온 지도

전국 만개일 예측 범위

선택 지역들의 만개일 변화 추이

GitHub 저장소

https://github.com/jungjae0/LEC-AgProgramming/tree/main/03_FloweringModels