Post

Earth Engine 데이터 다루기 with Python

Posted Updated
By Jung Jaeyoung 12 min read

본 문서는 Earth Engine Python API를 사용하여 ee.Image/Feature 등을 Pandas의 DataFrame 또는 Numpy의 Array 형태로 얻는 방법을 담고 있다.

Earth Engine과 Streamlit에서는 ImageCollection으로부터 직접 DataFrame을 생성하였다면, computeFeatures을 사용하여 바로 type을 바꾸는 방법을 담고 있다.

1. Library

필요한 library를 import 하고, earth engine 계정 인증과 초기화를 진행한다.

1
2
3
4
5
6
7
8

import altair as alt
import ee
import eerepr
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from mpl_toolkits.axes_grid1 import ImageGrid

1
2

# ee.Authenticate()
ee.Initialize()

2. Data

2-1. BASINS

WWF HydroSHEDS Basins Level 6

전세계의 500m 해상도의 수로 정보를 제공한다.

2-2. BOUNDARIES

FAO GAUL: Global Administrative Unit Layers 2015, First-Level Administrative Units

Global Administrative Unit Layers(GAUL, 전세계 행정구역 레이어)

2-3. CLIMATE

WorldClim Climatology V1

최저, 평균, 최고 기온과 강수량에 대한 월간 평균 지구 기후 데이터를 제공한다.

1
2
3

BASINS_ID = 'WWF/HydroSHEDS/v1/Basins/hybas_6'
BOUNDARIES_ID = 'FAO/GAUL/2015/level1'
CLIMATE_ID = 'WORLDCLIM/V1/MONTHLY'

BOUNDARIES_ID에서 대한민국을 필터링하여 대한민국 레이어(ko)를 얻고, ko를 사용하여 BASINS_ID에서 대한민국을 필터링한다.

1
2
3
4
5
6

basins = ee.FeatureCollection(BASINS_ID)
ko = ee.FeatureCollection(BOUNDARIES_ID).filter(
ee.Filter.eq('ADM0_NAME', 'Republic of Korea')
)

ko_basins = basins.filterBounds(ko)

CLIMATE_ID에서 강수량 밴드를 부분 집합화하고 개별 이미지를 단일 이미지로 쌓는다.

1
2
3
4
5
6
7
8
9

precip = ee.ImageCollection(CLIMATE_ID).select('prec')

months = precip.aggregate_array('month').getInfo()

band_names = [f'prec_month_{str(m).zfill(2)}' for m in months]

monthly_precip = ee.Image(precip.toBands().rename(band_names))

# monthly_precip.bandNames()

각 유역(ko_basins)의 평균 월 강수량을 계산한다.

1
2
3
4
5

ko_basins = monthly_precip.reduceRegions(
    collection=ko_basins,
    reducer=ee.Reducer.mean(),
    scale=1e3
)

3. FeatueCollection to Pandas’s DataFrame

ee.data.computeFeatures에서 'fileFormat': 'PANDAS_DATAFRAME'`을 사용하여 데이터를 pandas의 DataFrame 형태로 가져온다.

1
2
3
4

ko_basins_df = ee.data.computeFeatures({
    'expression': ko_basins,
    'fileFormat': 'PANDAS_DATAFRAME'
})

1
2

display(type(ko_basins_df))
ko_basins_df.head()

1

pandas.core.frame.DataFrame
geoCOASTDIST_MAINDIST_SINKENDOHYBAS_IDMAIN_BASNEXT_DOWNNEXT_SINKORDER...prec_month_03prec_month_04prec_month_05prec_month_06prec_month_07prec_month_08prec_month_09prec_month_10prec_month_11prec_month_12
0{'type': 'Polygon', 'coordinates': [[[127.4993...000040600039404060003940040600039401...58.02841499.99004593.973703162.829607263.427185232.806206143.81068452.77961244.20566924.069764
1{'type': 'Polygon', 'coordinates': [[[126.8624...000040600041004060004100040600041001...62.543293118.377280109.540524204.908318303.571153265.380813154.83321158.99214650.46605332.036306
2{'type': 'Polygon', 'coordinates': [[[126.4374...000040600045004060004500040600045001...61.486832105.254806101.688776178.860403259.184409229.985039148.62885953.99702852.83342733.797719
3{'type': 'Polygon', 'coordinates': [[[126.6916...000040600046804060004680040600046801...53.85280993.42782691.965696156.287028303.727885250.188197138.76981250.89650248.90372730.928487
4{'type': 'Polygon', 'coordinates': [[[127.4375...000040600050504060005050040600050501...51.46840188.61798496.304852150.415489338.487644286.903079168.18826756.98293449.86207129.780648

5 rows × 26 columns

3-1. melt

1
2
3
4
5
6
7

ko_basins_df = ko_basins_df.melt(
    id_vars=["HYBAS_ID"],
    value_vars=band_names,
    var_name="Month",
    value_name="Precipitation",
)
ko_basins_df.head()
HYBAS_IDMonthPrecipitation
04060003940prec_month_0127.502053
14060004100prec_month_0135.703600
24060004500prec_month_0136.438101
34060004680prec_month_0128.868651
44060005050prec_month_0128.169797

3-2. groupby & bar plot

groupby를 사용하여 각 유역의 평균 연간 강수량 그래프를 그린다.

1

ko_basins_df.groupby(['HYBAS_ID'])['Precipitation'].sum().plot.bar();

그림1

1
2
3
4
5
6

alt.Chart(ko_basins_df).mark_bar().encode(
    x=alt.X('HYBAS_ID:O'),
    y=alt.Y('Precipitation:Q', title='Precipitation (mm)'),
    color=alt.Color('Month', scale=alt.Scale(scheme='rainbow')),
    tooltip=alt.Tooltip(['HYBAS_ID', 'Precipitation', 'Month'])
).interactive()

4. FeatureCollection to GeoPandas’s GeoDataFrame

ee.data.computeFeatures에서 'fileFormat': 'GEOPANDAS_GEODATAFRAME'`을 사용하여 데이터를 geopandas의 GeoDataFrame 형태로 가져온다. 좌표계는 EPSG:5186로 지정해준다.

1
2
3
4
5
6
7
8

ko_basins_gdf = ee.data.computeFeatures({
    'expression': ko_basins,
    'fileFormat': 'GEOPANDAS_GEODATAFRAME'
})
ko_basins_gdf.crs = 'EPSG:5186'

display(type(ko_basins_gdf))
ko_basins_gdf.head()

1

geopandas.geodataframe.GeoDataFrame
geometryCOASTDIST_MAINDIST_SINKENDOHYBAS_IDMAIN_BASNEXT_DOWNNEXT_SINKORDER...prec_month_03prec_month_04prec_month_05prec_month_06prec_month_07prec_month_08prec_month_09prec_month_10prec_month_11prec_month_12
0POLYGON ((127.499 35.408, 127.504 35.402, 127....000040600039404060003940040600039401...58.02841499.99004593.973703162.829607263.427185232.806206143.81068452.77961244.20566924.069764
1POLYGON ((126.862 35.471, 126.864 35.467, 126....000040600041004060004100040600041001...62.543293118.377280109.540524204.908318303.571153265.380813154.83321158.99214650.46605332.036306
2POLYGON ((126.437 35.175, 126.438 35.168, 126....000040600045004060004500040600045001...61.486832105.254806101.688776178.860403259.184409229.985039148.62885953.99702852.83342733.797719
3POLYGON ((126.692 36.192, 126.694 36.191, 126....000040600046804060004680040600046801...53.85280993.42782691.965696156.287028303.727885250.188197138.76981250.89650248.90372730.928487
4POLYGON ((127.438 38.742, 127.437 38.710, 127....000040600050504060005050040600050501...51.46840188.61798496.304852150.415489338.487644286.903079168.18826756.98293449.86207129.780648

5 rows × 26 columns

1
2
3
4
5
6
7

ko_gdf = ee.data.computeFeatures({
    'expression': ko,
    'fileFormat': 'GEOPANDAS_GEODATAFRAME'
})

ko_gdf.crs = 'EPSG:5186'
ko_basins_gdf = ko_basins_gdf.clip(ko_gdf).to_crs(ko_gdf.crs)

4-1. prec_total

각 유역의 12개월 평균 강수량을 합산하여 prec_total열을 만든다.

1
2

ko_basins_gdf['prec_total'] = ko_basins_gdf[band_names].sum(axis=1)
ko_basins_gdf.head()
geometryCOASTDIST_MAINDIST_SINKENDOHYBAS_IDMAIN_BASNEXT_DOWNNEXT_SINKORDER...prec_month_04prec_month_05prec_month_06prec_month_07prec_month_08prec_month_09prec_month_10prec_month_11prec_month_12prec_total
6MULTIPOLYGON (((126.166 33.282, 126.165 33.283...100040600340104060034010040600340100...142.514948152.581449223.035944239.163145234.208038191.59062467.01404770.79292646.2378661594.097751
9MULTIPOLYGON (((125.209 34.685, 125.209 34.685...100040600041104060004110040600041100...125.082529126.071002215.532278254.175067242.603207157.44739054.74914950.12821028.3871791404.952484
2POLYGON ((126.438 35.179, 126.444 35.183, 126....000040600045004060004500040600045001...105.254806101.688776178.860403259.184409229.985039148.62885953.99702852.83342733.7977191307.807837
10MULTIPOLYGON (((126.010 35.326, 126.001 35.324...100040600045104060004510040600045100...95.28216991.067951152.975195258.768160220.711349137.52679152.84618454.32668436.8895451234.627188
11MULTIPOLYGON (((126.021 37.065, 126.014 37.063...100040600046904060004690040600046900...90.99958486.807500131.694396301.573913244.272398144.17022250.63988950.63540230.5377851234.594562

5 rows × 27 columns

4-2. 단계구분도

단계구분도를 만들고, 최소 강수량 지역은 빨간색 테두리로, 최대 강수량 지역은 분홍색 테두리로 표시해본다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# Define the choropleth map.
ax = ko_basins_gdf.plot(
    column='prec_total',
    cmap='viridis_r',
    vmin=ko_basins_gdf['prec_total'].min(),
    vmax=ko_basins_gdf['prec_total'].max(),
    legend=False,
    edgecolor='grey', linewidth=0.5
)

# Highlight the basin with the minimum annual precipitation: subset the geometry
# with the minimum precipitation total and then add it to the basin
# precipitation plot.
min_prec_gdf = ko_basins_gdf.loc[[ko_basins_gdf['prec_total'].idxmin()]]
min_prec_gdf.plot(ax=ax, color='none', edgecolor='red', linewidth=2)

max_prec_gdf = ko_basins_gdf.loc[[ko_basins_gdf['prec_total'].idxmax()]]
max_prec_gdf.plot(ax=ax, color='none', edgecolor='pink', linewidth=2)
# Add axis labels, a colorbar, and rotate x axis ticks.
ax.set_xlabel('Eastings [m]')
ax.set_ylabel('Northings [m]')
colorbar = plt.colorbar(ax.get_children()[0], fraction=0.03)
colorbar.set_label('Precipitation (mm)')
plt.xticks(rotation=45)

plt.show()

그림1

5. Image to Numpy Array

ee.data.computeFeatures에서 'fileFormat': 'NUMPY_NDARRAY'`을 사용하여 데이터를 Numpy 구조 배열로 계산된 이미지 형태로 가져온다.

1
2
3
4
5
6
7
8
9
10

monthly_precip_korea = monthly_precip.clipToBoundsAndScale(
    geometry=ko_basins, scale=1500
)

monthly_precip_npy = ee.data.computePixels({
    'expression': monthly_precip_korea,
    'fileFormat': 'NUMPY_NDARRAY'
})

monthly_precip_npy.shape

1

(609, 541)

이 작업은 다중 대역 이미지 데이터에 적합하다. 각 밴드가 열 이름인 배열 테이블로 생각할 수 있다. 또, 각 밴드가 서로 다른 데이터 유형을 가질 수 있다.

5-1. 한 장의 지도

1
2
3
4
5
6

names = monthly_precip_npy.dtype.names
print('field names:', names)

prec_month_10_arr = monthly_precip_npy['prec_month_10']
print('Selected array (band) shape:', prec_month_10_arr.shape)
plt.imshow(prec_month_10_arr, vmin=0, vmax=320);

1
2

field names: ('prec_month_01', 'prec_month_02', 'prec_month_03', 'prec_month_04', 'prec_month_05', 'prec_month_06', 'prec_month_07', 'prec_month_08', 'prec_month_09', 'prec_month_10', 'prec_month_11', 'prec_month_12')
Selected array (band) shape: (609, 541)

그림1

5-2. 시계열 강수량 지도

각 달의 평균 강수량을 그려 시계열로 표시해본다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# Set up the figure and grid.
fig = plt.figure(figsize=(20.0, 20.0))
grid = ImageGrid(
    fig,
    111,
    nrows_ncols=(4, 3),
    axes_pad=0.4,
    cbar_mode="single",
    cbar_location="right",
    cbar_pad=0.4,
    cbar_size="2%",
)

# Display each band to a grid cell.
for ax, name in zip(grid, names):
    ax.imshow(monthly_precip_npy[name], vmin=0, vmax=500)
    ax.set_title(name)

# Add colorbar.
colorbar = plt.colorbar(ax.get_children()[0], cax=grid[0].cax)
colorbar.set_label("Precipitation (mm)")

plt.show()

그림1

6. Soreted Data

6-1. FeatureCollection to DataFrame

속성 값으로 필터링하는 것이 가능하다. 한국의 유역으로 필터링하고, 하천 차수가 1이상인 유역만 DataFrame 또는 GeoDataFrame으로 가져온다. ORDER로 필터링을 하면 16개 중 6개만 남기 때문에 주석 처리하였다.

1
2
3
4
5
6
7
8

high_order_ko_basins_df = ee.data.listFeatures({
    'assetId': 'WWF/HydroSHEDS/v1/Basins/hybas_6',
    'region': ko.geometry().getInfo(),
    # 'filter': 'ORDER >= 1',
    'fileFormat': 'PANDAS_DATAFRAME'
})

high_order_ko_basins_df.head()
geoCOASTDIST_MAINDIST_SINKENDOHYBAS_IDMAIN_BASNEXT_DOWNNEXT_SINKORDERPFAF_IDSORTSUB_AREAUP_AREA
0{'type': 'Polygon', 'coordinates': [[[127.4993...00004060003940406000394004060003940142410231723296.423296.4
1{'type': 'Polygon', 'coordinates': [[[126.8624...0000406000410040600041000406000410014241043194833.04833.0
2{'type': 'Polygon', 'coordinates': [[[126.4374...0000406000450040600045000406000450014241063212860.92861.4
3{'type': 'Polygon', 'coordinates': [[[126.6916...0000406000468040600046800406000468014241083239776.89776.8
4{'type': 'Polygon', 'coordinates': [[[127.4375...00004060005050406000505004060005050142420032525245.325245.3
1
2
3
4
5
6
7
8
9

high_order_ko_basins_gdf = ee.data.listFeatures({
    'assetId': 'WWF/HydroSHEDS/v1/Basins/hybas_6',
    'region': ko.geometry().getInfo(),
    # 'filter': 'ORDER >= 1',
    'fileFormat': 'GEOPANDAS_GEODATAFRAME'
})

display(type(ko_basins_gdf))
high_order_ko_basins_gdf.head()

1

geopandas.geodataframe.GeoDataFrame
geometryCOASTDIST_MAINDIST_SINKENDOHYBAS_IDMAIN_BASNEXT_DOWNNEXT_SINKORDERPFAF_IDSORTSUB_AREAUP_AREA
0POLYGON ((127.49939 35.40779, 127.50356 35.402...00004060003940406000394004060003940142410231723296.423296.4
1POLYGON ((126.86250 35.47083, 126.86371 35.466...0000406000410040600041000406000410014241043194833.04833.0
2POLYGON ((126.43750 35.17500, 126.43824 35.167...0000406000450040600045000406000450014241063212860.92861.4
3POLYGON ((126.69167 36.19166, 126.69412 36.191...0000406000468040600046800406000468014241083239776.89776.8
4POLYGON ((127.43750 38.74167, 127.43715 38.710...00004060005050406000505004060005050142420032525245.325245.3
1
2
3
4
5
6
7
8
9
10
11

# Create an initial plot with the high river order watersheds.
ax = high_order_ko_basins_gdf.plot(edgecolor='purple', linewidth=1)

# Overlay the Washington state border for context.
ko_gdf.plot(ax=ax, color='none', edgecolor='black', linewidth=1)

# Set axis labels.
ax.set_xlabel('Eastings [degrees]')
ax.set_ylabel('Northings [degrees]')

plt.show()

그림1

6-2. Image to Numpy Array

ee.Image.clipToBoundsAndScale함수를 사용해 요청 영역과 스케일을 정의할 수 없다. 따라서 grid 매개변수를 사용해서 영역을 지정한다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

SCALE_FACTOR = 5

jan_mean_temp_npy = ee.data.getPixels({
    'assetId': 'WORLDCLIM/V1/MONTHLY/01',
    'fileFormat': 'NUMPY_NDARRAY',
    'grid': {
        'dimensions': {
            'width': 360 * SCALE_FACTOR,
            'height': 180 * SCALE_FACTOR
        },
        'affineTransform': {
            'scaleX': 1 / SCALE_FACTOR,
            'shearX': 0,
            'translateX': -180,
            'shearY': 0,
            'scaleY': -1 / SCALE_FACTOR,
            'translateY': 90
        },
        'crsCode': 'EPSG:4326',
    },
    'bandIds': ['tavg']
})

jan_mean_temp_npy.shape

1

(900, 1800)

배열에서 -9999보다 작은 값들은 Nan으로 대체한다. 배열의 모든 원소에 0.1을 곱해 값을 스케일링해 실제 온도로 변환한다.

1
2
3
4

jan_mean_temp_npy = jan_mean_temp_npy['tavg']

jan_mean_temp_npy = np.where(jan_mean_temp_npy < -9999, np.nan, jan_mean_temp_npy)
jan_mean_temp_npy = jan_mean_temp_npy * 0.1

1
2
3
4
5
6
7

fig = plt.figure(figsize=(10., 10.))
ax = plt.imshow(jan_mean_temp_npy, cmap='coolwarm', vmin=-40, vmax=40)

colorbar = plt.colorbar(ax, fraction=0.0235)
colorbar.set_label('Mean January Temp (°C)')

plt.show()

그림1

Earth Engine
Earth Engine
This post is licensed under CC BY 4.0 by the author.