기상 정보와 위성영상 기반의 식생지수를 사용한 머신러닝 기반 논 농작물 생산량 예측
프로젝트 개요
농작물 생산량 예측은 농업 자원 관리를 최적화하고, 시장 수급 조절 및 가격 예측을 통해 정책 수립에 중요한 역할을 한다. 농작물 생산량 예측을 위해서는 비료 사용, 토양 조성 등의 정보를 반영할 필요가 있으나, 개별 농가의 각종 데이터를 얻는 데는 한계가 있다. 이를 해결하기 위한 방안으로 위성영상 또는 무인 비행체 영상을 통해 실시간으로 넓은 지역의 식생을 파악할 수 있다. 또한 기상청을 통해 기상 정보를 얻을 수 있다. 본 프로젝트에서는 농작물 생산을 모니터링할 수 있는 식생정보와 농작물 생산에 영향을 미치는 요소 중 하나인 기상정보를 바탕을 농작물 생산량을 예측하고자 하였다.
재료 및 방법
1. 대상 지역, 기간, 농작물 선정과 생산량 데이터 취득
기상의 영향을 많이 받으며, 위성영상을 사용하여 식생을 모니터링할 수 있는 논 작물, 미곡, 두류, 잡곡, 맥류를 대상으로 하였다. 위성영상 데이터 취득에 많은 시간이 소요되고, 읍면동별로 생산량 자료를 제공하는 시군이 많지 않았다. 따라서 생산량 자료 제공 여부와 통계청의 ‘시군별 논밭별 경지면적’을 참고하여 전라남도 해남군, 전북특별자치도 김제시, 충청남도 당진시 3개 지역을 선정하였다. 또한 통계청에서 제공하는 생산량 데이터에서 세 지역의 공통 제공 기간인 2021년부터 2014년까지를 대상 기간으로 정하였다. 통계청 제공 읍면동별 연도별 미곡, 두류, 잡곡, 맥류 생산량 데이터를 취득하였다.
2. 식생지수 선정
식생지수는 식물의 상태, 분포 및 생장 상태를 평가하기 위해 사용되는 광학 원격 탐사 지수이다. 식생지수는 주로 위성/항공 영상을 분석하여 식생 상태를 정량적으로 평가하는 데 사용된다. 엽록소에 반사되는 특성을 가져 식생의 작황 추정에 적합한 NIR을 조합으로 하는 5개의 식생지수를 선정하였다.
| 지표 | 계산식 | 설명 |
|---|---|---|
| NDVI | (NIR - Red) / (NIR + Red) |
식생지수(Normalized Difference Vegetation Index), NIR(근적외선)과 Red(적색) 밴드를 이용해 식물의 건강 상태를 나타냄. |
| NDRE | (NIR - Red Edge) / (NIR + Red Edge) |
근적외선(NIR)과 빨간색 엣지(Red Edge) 밴드를 이용하여 식물의 엽록소 농도나 생육 상태를 평가. |
| GNDVI | (NIR - Green) / (NIR + Green) |
초록색 밴드(Green)와 근적외선(NIR) 밴드를 이용해 식물의 생장 상태를 평가. |
| RVI | NIR / Red |
반사율지수(Ratio Vegetation Index), NIR와 Red 밴드의 비율을 이용하여 식물의 생장 상태를 평가. |
| CVI | (NIR - Blue) / (Red + Blue) |
식물의 엽록소 농도나 건강 상태를 나타내는 지표로, NIR과 Blue, Red 밴드를 이용하여 계산. |
3. 위성영상 기반 식생지수 데이터 취득
Sentinel-2와 Landsat-8은 지구 관측을 위해 설계된 위성 시스템이다. 전 세계를 촬영하여 고해상도의 다중 스펙트럼 이미지를 제공한다. Sentinel-2는 MultiSpectralInstrument(MSI)를 탑재하고 있으며, 가시광선, 근적외선, 단파적외선 등 13개의 스펙트럼 밴드를 제공한다. Blue, Green, Red, Near Infrared 밴드는 10m 공간해상도, Red Edge는 20m 공간해상도이며, 5일을 주기로 촬영한다. Landsat8은 다양한 지표와 물체의 특성을 탐지하는 데 사용하는 OLI는 가시광선, 근적외선 등 9개의 스펙트럼 밴드를, 지표면의 열 방출을 측정하는 데 사용하는 Thermal Infrared Sensor(TIRS)는 2개의 열적외선 밴드를 제공한다. Panchromatic(흑백 이미지)는 15m, 가시광선, 근적외선, 단파적외선은 30m, 열적외선은 100m 공간해상도이며, 16일을 주기로 촬영한다. Sentinel-2는 2017년 3월 이후, Landsat8은 2013년 3월 이후부터 제공하기 때문에 2021~2019년은 Sentinel-2, 2018~2014년은 Landsat8를 사용했다. Google Earth Engine Python API를 사용하여 위성영상으로부터 식생지수 데이터를 취득하였다. 대한민국의 위성영상을 다룰 경우 ‘메모리 초과’ 에러가 발생한다. 따라서 읍면동 경계 json 데이터를 사용하여 읍면동 기준으로 개별 데이터를 취득하였다. 방법은 아래와 같다.
- 구름과 그림자 제거 구름과 구름 그림자를 제거하여 더 깨끗한 이미지를 얻었다. 구름 품질에 관한 정보를 담고 있는 밴드(Sentinel-2: QA60, Landsat-8: QA PIXEL)를 선택하고, 구름 마스크와 구름 그림자 마스크를 설정한 뒤, 구름과 그림자가 없는 부분을 식별한다. 구름 및 구름 그림자가 없는 부분을 결합하여 최종 마스크를 만든다. 이 마스크를 이미지에 적용하여 구름과 구름 이미지가 제거된 이미지를 반환한다. ```python qa = image.select(‘QA60’)
1.1. 구름과 권운 비트 마스크
cloud_bit_mask = 1 « 10 cirrus_bit_mask = 1 « 11
1.2. 마스킹 조건 정의
mask = (qa.bitwiseAnd(cloud_bit_mask).eq(0) .And(qa.bitwiseAnd(cirrus_bit_mask).eq(0)))
image.updateMask(mask)
2. 방사 보정
각 위성에서 획득한 원시 데이터는 디지털 수치(DN) 값으로 제공된다.
이러한 DN 값은 실제 물리적 의미를 갖지 않으므로, 이를 물리적 단위로 변환하는 작업이 필요하다.
센서 특성, 지형, 대기 효과 등에 의해 왜곡되는데,
이러한 오차를 보정계수를 이용해 수정하고 지형 지물에 대한 순수한 값을 구하는 작업이 필요하다.
```python
optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
image.addBands(optical_bands, None, True).addBands(thermal_bands, None, True)
- 식생지수 밴드 추가
NDVI, NDRE, GNDVI, RVI, CVI 밴드를 추가한다.
```python
1. 식생지수 산출
1.1. NDVI = (NIR - RED) / (NIR + RED)
ndvi = image.normalizedDifference([band_dct[‘NIR’], band_dct[‘RED’]]).rename(‘NDVI’)
1.2. NDRE = (NIR - RED_EDGE) / (NIR + RED_EDGE)
ndre = image.normalizedDifference([band_dct[‘NIR’], band_dct[‘RED_EDGE’]]).rename(‘NDRE’)
1.3. GNDVI = (NIR - GREEN) / (NIR + GREEN)
gndvi = image.normalizedDifference([band_dct[‘NIR’], band_dct[‘GREEN’]]).rename(‘GNDVI’)
1.4. RVI = NIR / RED
rvi = NIR.divide(RED).rename(‘RVI’)
1.5. CVI = (RED / GREEN^2) * NIR
cvi = image.expression( ‘(RED / GREEN ** 2) * NIR’, { ‘NIR’: image.select(band_dct[‘NIR’]), ‘GREEN’: image.select(band_dct[‘GREEN’]), ‘RED’: image.select(band_dct[‘RED’]) }).rename(‘CVI’)
2. 위성영상에 식생지수 밴드 추가
image.addBands([ndvi, cvi, ndre, gndvi, rvi])
4. 위성영상 재배열 및 저장
Sentinel-2와 Landsat8의 공간해상도가 다르다. 따라서 30m 해상도로 이미지로 재배열한다. 좌표 참조 시스템은 EPSG:4326으로 설정한다. 최종적으로 촬영일자별 각 식생지수 밴드에 대한 tif 파일을 저장한다.
```python
file_path = os.path.join(local_folder, f'Image_{date}.zip')
if not os.path.exists(file_path):
# 2.1. 해상도 = 30m, 좌표계 = WGS84
url = image.select(bands).getDownloadURL({
'scale': 30,
'crs': 'EPSG:4326',
'region': geometry
})
response = requests.get(url, stream=True)
with open(file_path, 'wb') as file:
for chunk in response.iter_content(chunk_size=8192):
file.write(chunk)
- 논 영역 필터링과 식생지수 데이터 취득 환경공간정보서비스에서 제공하는 2023년 중분류 토지피복도 shp 파일을 활용하여 논 영역을 필터링하였다. 각 지수별로 임의의 임계값을 설정하고, 평균값, 최솟값, 최댓값, 중앙값을 얻었다. 최종 산출물은 읍면동별 촬영일자별 식생지수를 담은 csv 파일이다. ```python fields = cover[cover[‘L2_NAME’] == ‘논’] if fields.crs != standard_crs: fields = fields.to_crs(standard_crs)
all_geometries = fields.geometry.unary_union geom_json = [all_geometries.geo_interface]
### 4. 기상 데이터 취득
1. 기상청 제공 ASOS
기상청 제공 ASOS API를 사용하여 일단위 기상 데이터를 취득하였다.
김제시의 경우 기상대가 없어 김제시와 가까운 부안군의 데이터를 취득하였다.
2. 월단위 기상 통계값
일별 평균 온도/습도는 웗별 평균 온도/습도, 일별 최고 온도는 월 최고 온도, 일별 최고 온도는 일 최저 온도, 일별 강수량은 총 강수량으로 처리하여 월단위 기상 통계 데이터를 취득하였다.
```python
df_each = df_each.groupby(['year', 'month']).agg({'tavg': 'mean',
'rainfall': 'sum',
'humid': 'mean',
'tmax': 'max',
'tmin': 'min'}).reset_index()
5. 최종 학습 데이터 생성
연도별 읍면동별 데이터를 얻고자 하였다. 위성영상의 경우 촬영주기는 일정하나, 연도별로 촬영일자가 다를 수 있기 때문에 월별로 번호를 매김으로써 월별 촬영 회차(3월 1회차, 3월 2회차, …)로 수정하였다. 최종적인 데이터는 아래와 같으며, csv 파일이다.
df['date'] = df['date'].astype(str).str.replace('.zip', '')
df['year'] = df['date'].str[:4]
df['month'] = df['date'].str[4:6]
df = df[(df['month'].astype(int) > 2) & (df['month'].astype(int) < 11)]
df['yearmonth'] = df['year'] + df['month']
df['day'] = df['date'].str[6:]
df['seq'] = df.groupby(['yearmonth', '행정구역']).cumcount() + 1
df['seq'] = df['seq'].apply(lambda n: chr(64 + n))
df['seq'] = df['month'] + df['seq']
df = df.drop(columns=['yearmonth', 'date', 'month', 'day'])
6. 더미변수 추가와 데이터스케일링
농작물 분류를 더미변수로 추가하여 해당 정보를 고려할 수 있도록 하였다. Standard Scaler를 사용하여 평균이 0이고, 분산이 1이 데이터를 표준화하여 데이터의 분포를 중앙에 맞추고, 서로 다른 특성 간의 스케일을 동일하게 하였다. 스케일링을 통해 모든 특성이 균등하게 고려되고, 안정적인 학습이 될 수 있도록 하였다.
7. 머신러닝 모델과 모델 평가
정교한 부스팅 방법과 뛰어난 성능을 보이는 XGBoost, 노이즈에 강하며 쉽고 직관적인 RandomForest 모델을 선정하였다. 미곡, 맥류, 두류, 잡곡 전체 세부 작물에 대한 2021년 이전의 생산량(M/T, ‘단위생산량’이 아닌 ‘생산량’)으로 2021년 생산량을 예측하였다. 평가 지표로는 R²와 RMSE를 사용하였다. Python scikit-learn 패키지를 사용하여 모델 학습과 평가를 진행하였다.
결과
1. 구름 마스킹 결과
위성영상에서 구름이 덮인 부분은 지표면의 실제 정보를 가리기 때문에, 이 부분을 그대로 사용하면 식생을 제대로 반영한다고 볼 수 없다.
구름은 빛을 산란, 반사하기 때문에 구름이 많은 영역은 다른 영역과 비교해 명도가 크게 다를 수 있으며, 이는 일관성을 확보하는 데 방해요소가 된다.
따라서 위성영상에서 구름을 처리하는 과정이 필요하다.
익산에 위치한 논 A를 대상으로 각 위성별로 구름 마스킹 여부에 따른 밀 생육단계별 식생지수 변화를 비교한 결과이다.
밀의 생장과정에서 점진적으로 녹색 면적이 증가하다가, 알곡이 형성되는 과정을 거치면서 녹색 면적이 줄어들기 때문에 식생지수 변화 그래프에서 언덕과 같은 모습을 보이는 것이 이상적이다.
두 위성 모두 구름 마스킹 전(점선)에는 큰 변동을 보인다.
Sentinel-2에서는 구름 마스킹에 따라 크게 개선되는 효과가 없었지만, Landsat-8의 경우 구름 마스킹을 함(실선)으로써 큰 변동을 보이지 않고 기대하는 경향을 따른다.
2. 연도별 기상 변화
세 지역 모두 7월과 8월에 강수량이 집중되는 경향을 보인다. 이는 한국의 여름철 장마와 태풍의 영향으로 보인다.
특히 2014년, 2017년, 2021년의 7월과 8월 강수량이 다른 연도에 비해 월등히 높다. 6월과 9월에도 약간의 강수량 증가를 보이지만, 나머지 월들은 강수량이 적고 일정한 패턴을 보이지 않는다.
세 지역 모두 월별 기온 변화 패턴이 유사하다. 모든 지역에서 1월과 2월에 기온이 가장 낮고, 7월과 8월에 가장 높다.
연도별로 큰 차이는 없으나, 해를 거듭할수록 기온 상승이 나타나는 경향이 있다. 전반적으로 기온 상승 추세를 볼 때, 기후 변화의 영향을 받았을 가능성이 있다.
3. 연도별 농작물 생산량 변화
미곡(노란색)은 세 지역에서 가장 많이 생산되는 작물이다.
해남군에서 변동폭이 가장 크고, 김제시와 부안군에서 비교적 일정한 생산량을 유지하고 있다.
맥류(빨간색)는 세 지역 모두에서 점진적인 증가세를 보이나 부안군을 제외한 지역에서는 2021년 감소하였다.
잡곡(파란색)은 세 지역 모두에서 생산량이 매우 적고 일정하게 유지되고 있다.
두류(녹색)는 세 지역 모두에서 소폭 증가하는 경향을 보이며, 김제시에서의 생산량 증가폭이 상대적으로 컸으나, 2021년 감소하였다.
4. 농작물 생산량 예측 결과
RandomForest 모델은 실제 값을 73.06% 설명할 수 있으며, RMSE는 1267.2878 M/T이다.
XGBoost 모델은 실제 값을 77.28% 설명할 수 있으며, RMSE는 1163.7392 M/T이다.
두 모델 모두 예측 값과 실제 값이 대체로 직선에 가까운 형태를 이루고 있다.
훈련 데이터(파란색 점)와 테스트 데이터(노란색 점)가 비교적 고르게 분포되어 있으나, 두 모델 모두 고값에서 예측 오차가 크게 나타나고 있다.
전반적으로 XGBoost 모델이 RandomForest 모델보다 조금 더 나은 예측 성능을 보인다.
그러나, 추가적인 모델 튜닝이나 다른 변수 고려, 추가 데이터 확보 등을 통해 모델을 보완할 필요가 있다.