solar radiation - time and space

2021-11-09

2025-12-30

geopandas, matplotlib, proj, python

한국에너지기술연구원은 천리안 위성 기반 전국 전일사량 데이터를 제공합니다.
천리안 1호 위성 데이터는 2012년 1월부터 2019년 12월까지,
천리안 2호 위성 데이터는 2019년 9월부터 2020년 8월까지 데이터를 제공합니다.
이 데이터 중 천리안 1호 데이터를 받아 분석합니다.

1. 데이터 다운로드

공공데이터포털 신재생자원지도 데이터: 태양자원 천리안 1호 수평면 전일사량

다운로드를 클릭하면 189 MB짜리 .csv 파일이 다운로드 됩니다.
파일 이름은 한국에너지기술연구원_신재생자원지도데이터_태양자원_천리안1호_수평면전일사량_20191231.csv입니다.

2. 데이터 분석

2.1. 기본 설정

주피터 노트북을 열고 기본 설정을 합니다.

%matplotlib inline

import matplotlib.pyplot as plt
import seaborn as sns
import geopandas as gpd
import numpy as np
import pandas as pd

# 시각화 설정
sns.set_context("talk")
sns.set_style("white")

# 한글 사용 설정
plt.rcParams['font.family']=['NanumGothic', 'sans-serif']
plt.rcParams['axes.unicode_minus'] = False

2.2. 데이터 불러오기

작업폴더 아래 data폴더를 만들고 다운받은 파일을 넣습니다.
pandas를 이용해 데이터를 불러옵니다.

df = pd.read_csv("./data/한국에너지기술연구원_신재생자원지도데이터_태양자원_천리안1호_수평면전일사량_20191231.csv", 
                 encoding="euc-kr")  # 한글 인코딩
df.to_pickle("./data/KIER.pkl")      # 백업
df.head()

위도와 경도, 그리고 시간별 숫자가 있습니다.
위도와 경도는 지구상의 지점을 지칭하고 시간별 숫자는 월별 평균 일사량을 구합니다.
직진에 의해 태양으로부터 전달되는 일사량과 구름에 의해 산란된 값을 합친 값입니다.
산란된 일사량을 산란일사라고 하는데, 에너지기술연구원에서는 도메인 지식 기반으로 계산하여 추론합니다.
데이터 형상을 확인합니다.

df.shape

실행 결과

1	(291924, 98)

약 100개 가량의 열과 30만개 행으로 구성되어 있습니다.
98개 행은 12개월 x 8년 + 위도, 경도입니다.

3. 시간에 대한 변화

3.1. 광주광역시 데이터 추출

광주광역시의 좌표는 북위 35.1595º, 동경 126.8526º 입니다.
광주광역시와 가장 가까이 위치한 지점의 좌표를 탐색합니다.
위도와 경도 모두 가까워야 하므로 제곱의 합이 최소인 지점을 구합니다.

1
2
3

# 광주광역시 위경도 : 35.1595° N, 126.8526° E
idx = np.argmin(np.power(df["lon"]-35.1595, 2) + np.power(df["lat"]-126.8526, 2))
df.iloc[idx]

시간에 대한 일사량을 확인하려면 시간과 일사량 데이터가 필요합니다.
월과 일사량 데이터를 가져옵니다.

1 2	months = list(df.columns)[2:] irrs = df.iloc[idx][2:]

추출한 데이터를 그림으로 그립니다.
ax.set_title()명령으로 데이터 이름을 제목으로 붙입니다.

fig, ax = plt.subplots(figsize=(10, 4))

# 변화율 시각화
ax.plot(months, irrs)

# 제목 추가
ax.set_title("수평면 전일사량 [kWh/(m$^2 \cdot$day)]", pad=10)

3.2. 범위 및 눈금 간격 지정

x축 눈금이 빼곡하게 적혀서 보이지 않습니다.
MultipleLocator를 사용해서 1월 데이터에만 눈금을 답니다.
ax.set_xlim()을 사용해서 x축 범위를 지정합니다.

from matplotlib.ticker import MultipleLocator

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(months, irrs)
ax.set_title("수평면 전일사량 [kWh/(m$^2 \cdot$day)]", pad=10)

# x축 범위 지정
ax.set_xlim(0, len(months)-1)

# x축 눈금 간격 지정
ax.xaxis.set_major_locator(MultipleLocator(12))

3.3. 눈금 글자 크기 축소

글자를 조금 줄이면 좋을 것 같습니다.
일사량 데이터이므로 따뜻한 느낌을 주는 오렌지색으로 바꿉니다.
1월마다 눈금을 만들어 올립니다.

from matplotlib.ticker import MultipleLocator

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(months, irrs, c="orange")
ax.set_title("수평면 전일사량 [kWh/(m$^2 \cdot$day)]", pad=10)
ax.set_xlim(0, len(months)-1)
ax.xaxis.set_major_locator(MultipleLocator(12))

# x축과 y축 눈금 글자 크기 축소
ax.tick_params(axis="both", labelsize="small")

# x축 눈금에서 grid 생성
ax.grid(axis="x")

3.4. 선 아래 채색

조금 심심하다 싶습니다.
ax.stackplot()명령을 사용해 선 아래를 색으로 칠합니다.

from matplotlib.ticker import MultipleLocator

fig, ax = plt.subplots(figsize=(12, 4))

# stackplot
ax.stackplot(months, irrs, colors="orange", ec="brown", alpha=0.5)

ax.set_title("수평면 전일사량 [kWh/(m$^2 \cdot$day)]", pad=10)
ax.set_xlim(0, len(months)-1)
ax.xaxis.set_major_locator(MultipleLocator(12))
ax.tick_params(axis="both", labelsize="small")
ax.grid(axis="x")

4. 공간 분포

4.1. 데이터 정돈

한반도를 포함한 주변 영역의 데이터를 2차원 공간에 그립시다.
데이터프레임 컬럼 이름이 너무 깁니다.
단위는 따로 적어두면 됩니다. 불필요한 부분을 제거합시다.
위도와 경도도 영어로 바꿉니다.

cols = list(df.columns)
cols_short = ["lon", "lat"] + [c.split(" ")[0] for c in cols[2:]]
df.columns = cols_short
df.head()

4.2. scatter plot

2차원 공간에 데이터를 점으로 찍어서 그림을 그립니다.
2012년 1월 데이터만 그립니다.

1 2	fig, ax = plt.subplots() ax.scatter(df["lat"], df["lon"], c=df["2012-01"], cmap="inferno")

Matplotlib에서 제공하는 그림 크기가 가로가 길고 세로가 작아 짓눌린 모양입니다.
세로로 긴 도화지를 깔고 종횡비(aspect ratio)를 1:1로 맞춥니다.
실제 형상을 반영하기 위해서입니다.

1
2
3

fig, ax = plt.subplots(figsize=(5, 6.67), constrained_layout=True)
ax.scatter(df["lat"], df["lon"], c=df["2012-01"], cmap="inferno")
ax.set_aspect("equal")

4.3. 좌표계 변환

위쪽이 넓고 아래가 좁은 모양입니다.
그림 아래와 왼쪽에 놓인 숫자는 경도와 위도입니다.
경도와 위도는 곡면상에 놓인 그림이기 때문에 이렇게 그리면 왜곡이 발생합니다.
pyproj 라이브러리를 사용해 UTM-K 좌표계로 변경합니다.

from pyproj import Transformer

transformer = Transformer.from_crs("epsg:4326", "epsg:5178")
coord_UTMK = transformer.transform(df["lon"], df["lat"])
coord_UTMK

변환된 좌표계를 기존 데이터에 붙여 데이터프레임을 만듭니다.

1
2
3

df_u = pd.DataFrame(data=np.array(coord_UTMK).T, columns=["y", "x"])
df_u = pd.concat([df_u, df.iloc[:, 2:]], axis=1)
df_u.head()

4.4. UTM-K 좌표계 scatter plot

변환된 좌표계를 사용해 scatter plot을 그립니다.

1
2
3

fig, ax = plt.subplots(figsize=(5, 6.67), constrained_layout=True)
ax.scatter(df_u["x"], df_u["y"], c=df_u["2012-01"], cmap="inferno")
ax.set_aspect("equal")

이제 그림이 정상적으로 보입니다.
경기 동북부~강원도 서부 지역을 확대해 봅시다.

fig, ax = plt.subplots(figsize=(5, 6.67), constrained_layout=True)
ax.scatter(df_u["x"], df_u["y"], c=df_u["2012-01"], cmap="inferno")
ax.set_aspect("equal")
ax.set_xlim(1e6, 1e6+10000)
ax.set_ylim(2e6, 2e6+10000)

확대해보니 점만 가득합니다.
scatter plot이라 점만 찍어 그린 것입니다.
위쪽 그림을 보니 아래쪽에 굵은 펜이 지나간 듯한 흔적이 보이기도 합니다.
조금 더 넓게 봅니다.

fig, ax = plt.subplots(figsize=(5, 6.7), constrained_layout=True)
ax.scatter(df_u["x"].iloc[:1000], 
           df_u["y"].iloc[:1000], 
           c=df_u["2012-01"].iloc[:1000], cmap="inferno")

가로 세로도 맞지 않습니다.
원래 데이터도 좀 이랬고, 좌표 변환의 여파이기도 합니다.

4.5. 2D image 변환

이처럼 x와 y가 쌍으로 얽힌 데이터를 2D로 표현하는 가장 정석적인 방법은 meshgrid를 만드는 것입니다.
그러나 이는 x, y 데이터가 격자모양으로 예쁘게 정리되었을 때나 가능합니다.
위성 영상은 x, y 위치가 조금씩 틀어져있기 때문에 grid를 만들기 적절치 않습니다.
이럴 때 triplot, tricontourf, tricontour를 사용하기 좋습니다.
명령어의 기능을 확인하기 위해 30만개 데이터 중 100개만 추출합니다.
그리고 같은 데이터를 세 가지 명령을 사용해 그립니다.

df_us = df_u.sample(100)

# visualize
fig, axs = plt.subplots(ncols=3, figsize=(9, 6.67), sharex=True, sharey=True, constrained_layout=True)
axs[0].triplot(df_us["x"], df_us["y"], lw=0.5)
axs[1].triplot(df_us["x"], df_us["y"], lw=0.5, c="w")
axs[2].triplot(df_us["x"], df_us["y"], lw=0.5, c="k", alpha=0.5)

axs[1].tricontourf(df_us["x"], df_us["y"], df_us["2012-01"], cmap="inferno", levels=5)
axs[2].tricontour(df_us["x"], df_us["y"], df_us["2012-01"], cmap="inferno", levels=5)

for ax, title in zip(axs, ["triplot", "tricontourf", "tricontour"]):
    ax.set_aspect("equal")
    ax.set_title(title, pad=20)

ax.triplot()은 주어진 (x, y) 데이터끼리 연결해 삼각형 mesh를 만듭니다.
ax.tricontourf()는 (x, y) 위치데이터에 node의 값 데이터를 적용해 삼각형 내 interpolation을 합니다.
levels=5로 지정했기 때문에 공간이 다섯 개의 색상 구간으로 나뉘었고 색상 사이 경계선이 보입니다.
ax.tricontour()는 이 경계선에 contour line을 그립니다.
이 방식을 사용해 2D에 분포된 데이터를 그림으로 표현할 수 있습니다.
일사량 데이터를 보이는데 ax.triplot()은 의미가 없으니 나머지 둘만 그립니다.
level도 10개로 설정해 더 많은 구간을 표현합니다.

fig, axs = plt.subplots(ncols=2, figsize=(6, 6.67), sharex=True, sharey=True, constrained_layout=True)

axs[0].tricontourf(df_u["x"], df_u["y"], df_u["2012-01"], cmap="inferno", levels=10)
axs[1].tricontour(df_u["x"], df_u["y"], df_u["2012-01"], cmap="inferno", levels=10)

for ax, title in zip(axs, ["tricontourf", "tricontour"]):
    ax.set_aspect("equal")
    ax.set_title(title, pad=20)

4.6. 지도 윤곽선 표시 (1) geopandas dataset

geopandas 라이브러리에는 세계 지도와 주요 도시가 들어있습니다.
이를 활용해 한반도와 주변국의 윤곽선을 그립니다.

import geopandas as gpd

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

# 국경 데이터
world.head(3)

1 2	# 도시 데이터 cities.head(3)

공간을 마련하고 한반도와 서울을 표시합니다.
geopandas 자체 기능을 사용해 국경을 그립니다.

1
2
3

fig, ax = plt.subplots()
world[world["name"] == "South Korea"].plot(ax=ax)
cities[cities["name"]=="Seoul"].plot(ax=ax, c="yellow")

자, 그런데 x와 y축 눈금을 보면 위도와 경도입니다.
앞에서 우리는 UTM-K 좌표계로 데이터 전체를 변환했습니다.
역시 geopandas 자체 기능을 사용해 데이터를 변환하고 다시 그립니다.

# 좌표계 변환
world = world.to_crs("EPSG:5178")
cities = cities.to_crs("EPSG:5178")

# 시각화
fig, ax = plt.subplots()
world[world["name"] == "South Korea"].plot(ax=ax)
cities[cities["name"]=="Seoul"].plot(ax=ax, c="yellow")

간편한 명령으로 손쉽게 좌표계가 변환되었습니다.
일사량 데이터와 겹쳐봅니다.

fig, ax = plt.subplots(figsize=(5, 6.67), constrained_layout=True)

ax.tricontourf(df_u["x"], df_u["y"], df_u["2012-01"], cmap="inferno", levels=10)
ax.set_aspect("equal")

world[world["name"] == "South Korea"].plot(ax=ax, fc="none", ec="w")
world[world["name"] == "North Korea"].plot(ax=ax, fc="none", ec="w")
world[world["name"] == "Japan"].plot(ax=ax, fc="none", ec="w")

ax.set_xlim(df_u["x"].min(), df_u["x"].max())
ax.set_ylim(df_u["y"].min(), df_u["y"].max())

디테일이 영 아쉽습니다.
심지어 제주도도 없는데 혹시나 싶어 북한과 일본까지 다 끌어와도 그려지지 않습니다.

4.7. 지도 윤곽선 표시 (2) shape file

외부 저장소에서 시도 경계선 데이터를 가져옵니다.
여기에서 다운로드받을 수 있습니다.
geopandas.read_file()명령으로 열어 확인합니다.

1 2	df_shp = gpd.read_file("./CTPRVN_202101/TL_SCCO_CTPRVN.shp", encoding="euc-kr") df_shp.head()

이번엔 그림으로 확인합니다.
특별시, 광역시, 도의 경계선이 선명하게 드러납니다.

1 2	ax = df_shp.plot() fig = ax.figure

4.8. 일사량 데이터와 overlap

아까 살펴본 2012년 1월 데이터에 겹쳐봅니다.

fig, ax = plt.subplots(figsize=(5, 6.67), constrained_layout=True)

ax.tricontourf(df_u["x"], df_u["y"], df_u["2012-01"], cmap="inferno", levels=10)
ax.set_aspect("equal")

df_shp.plot(ax=ax, fc="none", ec="w", lw=1, alpha=0.5)

ax.set_xlim(df_u["x"].min(), df_u["x"].max())
ax.set_ylim(df_u["y"].min(), df_u["y"].max())

그럴싸합니다!
추운 겨울을 그려봤으니 6개월 지난 7월 그림을 그려봅니다.

fig, ax = plt.subplots(figsize=(5, 6.67), constrained_layout=True)

ax.tricontourf(df_u["x"], df_u["y"], df_u["2012-07"], cmap="inferno", levels=10)
ax.set_aspect("equal")

df_shp.plot(ax=ax, fc="none", ec="w", lw=1, alpha=0.5)

ax.set_xlim(df_u["x"].min(), df_u["x"].max())
ax.set_ylim(df_u["y"].min(), df_u["y"].max())

데이터가 바뀌어 색상 분포가 달라졌습니다.
그런데 뭔가 수상합니다.
1월 데이터는 아주 밝았는데 7월 데이터는 어두침침합니다.
여름 일사량이 겨울보다 적다는 게 정상일까요?
확인을 위해 각각 colorbar를 붙여봅니다.

fig, axs = plt.subplots(ncols=2, figsize=(10, 6.67), constrained_layout=True, sharex=True, sharey=True)

im0 = axs[0].tricontourf(df_u["x"], df_u["y"], df_u["2012-01"], cmap="inferno", levels=10)
im1 = axs[1].tricontourf(df_u["x"], df_u["y"], df_u["2012-07"], cmap="inferno", levels=10)

for ax, im, month in zip(axs, [im0, im1], ["2012-01", "2012-07"]):
    df_shp.plot(ax=ax, fc="none", ec="w", lw=1, alpha=0.5, zorder=2)
    im = ax.tricontourf(df_u["x"], df_u["y"], df_u[month], cmap="inferno", levels=10)
    ax.set_aspect("equal")
    ax.set_title(month, pad=20)
    ax.set_xlim(df_u["x"].min(), df_u["x"].max())
    ax.set_ylim(df_u["y"].min(), df_u["y"].max())

    plt.colorbar(im, ax=ax)

2월은 최대값이 3을 넘지 않는데 비해 7월은 5.4에 달합니다.
최소값도 0.6과 1.8로 큰 차이가 납니다.
두 그림이 다른 기준으로 그려져서 문제가 생긴 것입니다.
위 코드의 6번째 줄에 colormap의 최소값vmin과 최대값vmax를 똑같이 지정해서 같은 기준으로 색을 입힙니다.

1	im = ax.tricontourf(df_u["x"], df_u["y"], df_u[month], cmap="inferno", levels=10, vmin=0, vmax=6)

확실히 1월에 비해 7월의 햇살이 더 쨍쨍합니다.

5. 시간과 공간 데이터 시각화

공간 데이터는 1개월마다 평균치가 저장되므로 이를 늘어놓으면 시간에 따른 변화를 볼 수 있습니다.
2012년 1월부터 12월까지 공간데이터를 for 문을 사용해 반복합니다.
공통으로 적용된 colorbar는 12개 Axes의 우측에 ScalarMappable()을 사용해 놓습니다.

from matplotlib import cm

fig, axes = plt.subplots(ncols=6, nrows=2, figsize=(10, 6), constrained_layout=True, sharex=True, sharey=True)

axs = axes.ravel()
for i, ax in enumerate(axs, 1):
    month = f"2012-{i:02d}"
    df_shp.plot(ax=ax, fc="none", ec="w", lw=1, alpha=0.5, zorder=2)
    im = ax.tricontourf(df_u["x"], df_u["y"], df_u[month], cmap="inferno", levels=10, vmin=0, vmax=6)
    ax.set_aspect("equal")
    ax.set_title(month, pad=10)
    ax.set_xlim(df_u["x"].min(), df_u["x"].max())
    ax.set_ylim(df_u["y"].min(), df_u["y"].max())
    ax.axis(False)

# colorbar
cbar = cm.ScalarMappable(cmap="inferno")
cbar.set_clim(0, 6)
plt.colorbar(cbar, ax=axes[:,-1], fraction=0.2, pad=0.15)

# figure title
fig.suptitle("수평면 전일사량 [kWh/(m$^2 \cdot$day)]", color="gray", fontsize="medium", fontweight="bold")

도움이 되셨나요? 카페인을 투입하시면 다음 포스팅으로 변환됩니다

PythonData Science

1. 데이터 다운로드

2. 데이터 분석

2.1. 기본 설정

2.2. 데이터 불러오기

3. 시간에 대한 변화

3.1. 광주광역시 데이터 추출

3.2. 범위 및 눈금 간격 지정

3.3. 눈금 글자 크기 축소

3.4. 선 아래 채색

4. 공간 분포

4.1. 데이터 정돈

4.2. scatter plot

4.3. 좌표계 변환

4.4. UTM-K 좌표계 scatter plot

4.5. 2D image 변환

4.6. 지도 윤곽선 표시 (1) geopandas dataset

4.7. 지도 윤곽선 표시 (2) shape file

4.8. 일사량 데이터와 overlap

5. 시간과 공간 데이터 시각화