본문으로 건너뛰기

위성이미지 처음 다지기: 로그인 없이 산 사진 받아 식생 지수까지 계산하는 7단계

· 약 16분
Datapopcorn CEO / AI automation educator

위성데이터는 어렵지 않습니다. 어려운 것은 위성이미지가 숫자가 아니라 그림이라는 착각입니다. 이 글은 로그인 없이 위성사진 한 장을 내 컴퓨터로 받아서, 거기서 숫자를 뽑는 방법까지 7단계로 정리한 안내서입니다.

이 글에서 사용하는 모든 명령과 모든 수치는 2026년 10월 3일에 직접 실행한 결과입니다. 되돌아온 오류 메시지도 그대로 옮겼습니다.

이 글은 이런 분께 필요합니다​

  • 필요한 배경: Python에서 함수를 하나 호출하고 리스트를 다루는 정도. 변수, 반복문, 함수 정의만 알면 충분합니다.
  • 필요한 도구: Python 3, 터미널, 인터넷. 계정은 하나도 필요 없습니다.
  • 필요한 시간: 처음이라면 20~30분. 7단계를 그대로 따라 하면 끝에 result.png가 생깁니다.
  • 이 글에서 하지 않는 것: 위성사진을 머신러닝으로 분류하거나, 학습 모델을 만드는 일은 이 글의 범위 밖입니다. 위 글은 "사진을 받고, 색을 만들고, 지수를 계산하고, 비교한다"까지만 다룹니다.

이 글에서 쓰는 프로그램​

글에서 말하는 캡처는 모두 2026년 10월 3일에 실제로 실행한 화면이나 그 결과물을 그대로 옮긴 것입니다.

단계쓰는 프로그램하는 일
Step 1터미널 (macOS 기본 Terminal)pip install 로 패키지 설치
Step 2~7VS Code + Python 3위성이미지 받기, 색 만들기, 지수 계산, 비교
Step 3 확인용브라우저 (Chrome)주소창에 URL 을 넣어 API 응답을 눈으로 확인
그림 확인용미리보기 (macOS Preview)내려받은 PNG 파일을 열어 눈으로 보기

추가로 쓰는 라이브러리는 셋입니다. pillow는 이미지를 열고 숫자 배열로 바꾸고, numpy는 화소 단위 계산을 하고, matplotlib은 지수를 색으로 칠해서 그림을 그립니다.

전체 구조를 먼저 보면​

위성데이터 작업은 크게 네 층위입니다. 처음에는 이 구조가 보이지 않아서 어디서 막힐지 모릅니다.

층위하는 일이 글의 단계
찾기어떤 위성사진이 내 관심 지역을 덮고 있는가Step 3
받기그중 한 장을 내 컴퓨터로 내려받는다Step 4
만들기밴드를 조합해 사람이 보는 색과 지수를 만든다Step 5, 6
비교하기두 시점의 사진을 겹쳐 변화를 찾는다Step 7

Step 1·2는 준비와 좌표 계산입니다. 자꾸 좌표가 나오는데, 위성데이터는 사진이 격자 위에만 존재하기 때문입니다.

Step 1. 터미널에서 준비물을 설치합니다​

프로그램: 터미널

위성사진은 보통 tif나 jp2 같은 전용 형식이라 예제 파일 뷰어로는 열지 않습니다. 이 글에서는 PNG 타일로 렌더링받는 방식을 씁니다. Pillow가 열 수 있어서 확인하기 편합니다.

터미널을 열고 아래 한 줄을 넣고 엔터를 치면 끝입니다.

pip install pillow matplotlib numpy

오류 없이 세 개가 설치됩니다.

Step 2. 볼 위치를 좌표로 정하고, 타일 번호로 바꿉니다​

프로그램: VS Code + Python

위성사진은 지구를 정사각형 격자로 자른 조각(타일)로 관리합니다. 그래서 "서울시청 앞" 같은 지명은 줄곧 타일 번호로 바꿔야 합니다.

웹 mercator 규격을 쓰면 변환은 수식 두 줄이면 됩니다. z는 확대 수준이고, z가 1 올라갈 때마다 타일 수가 네 배가 됩니다.

import math

def lonlat_to_tile(lon, lat, z):
n = 2 ** z
x = int((lon + 180) / 360 * n)
r = math.radians(lat)
y = int((1 - math.log(math.tan(r) + 1 / math.cos(r)) / math.pi) / 2 * n)
return x, y

# 설악산 대관령
LON, LAT = 128.551, 37.952
x, y = lonlat_to_tile(LON, LAT, 13)
print(x, y) # 7021 3161

VS Code에서 이 파일을 만들고 돌리면 터미널에 이렇게 나옵니다.

7021 3161

설악산 대관령이 떨어지는 z13 타일 위치. 왼쪽은 격자 위 좌표, 오른쪽은 그 타일이 실제로 담는 범위 요약입니다.

z를 13으로 두면 타일 한 장은 약 3.8km 정사각이 되고, 화소 크기는 약 15m입니다. 관심 지역을 훑어보려면 z 13~14가 무난합니다.

여기서 초보자가 제일 많이 헷갈리는 지점을 짚어 둡니다. Step 3에서 52SDH라는 MGRS 번호를 보게 되는데, 그건 약 100km짜리 타일입니다. 여기서 만든 7021, 3161은 그 안에서 잘라낸 약 3.8km짜리 웹 타일입니다. 두 개는 이름만 비슷하고 크기가 100배 넘게 다릅니다. 이 구분을 모르면 Step 4에서 사진 크기가 이상하다고 생각하게 됩니다.

z를 얼마나 크게 하느냐가 곧 해상도입니다. Sentinel-2는 원본 화소가 10m지만, 이 타일 경로에서는 타일 크기 자체가 최대 256×256 화소라 z가 커질수록 화소가 늘어나고 뭉개집니다.

Step 3. 그 위치를 덮는 위성사진을 찾습니다​

프로그램: Python (검색) + 브라우저 (눈으로 확인)

여기가 이 글에서 가장 많이 막히는 지점입니다. 검색을 아무리 잘 해도, 검색 결과가 내 좌표를 덮는다는 보장은 없습니다.

STAC 검색은 좌표 범위와 시간이 겹치는 장면을 전부 돌려줍니다. 그중 일부는 내 지점이 타일 경계 바깥이어서, 사진을 받으면 아무 것도 없는 검은 화면이 나옵니다. 이건 실측에서 실제로 겪은 상황입니다.

설악산 대관령(37.952N)을 검색했을 때 2026년 5월 장면이 두 개 나왔는데, 그중 구름이 가장 적은 장면의 북쪽 경계가 37.947도였습니다. 대관령은 37.952도라 타일 밖이었습니다. 구름 0.004%인 거의 완벽한 장면인데, 볼 수 있는 영역이 하나도 없는 상태였습니다.

그래서 검색 결과마다 직접 그 위치를 덮는지 확인하는 코드가 반드시 필요합니다.

import json, urllib.request

def find_scene(date_str, lon, lat, max_cloud=20):
day = f"{date_str}T00:00:00Z/{date_str}T23:59:59Z"
body = json.dumps({
"collections": ["sentinel-2-l2a"],
"bbox": [lon - 0.1, lat - 0.05, lon + 0.1, lat + 0.05],
"datetime": day,
"limit": 60,
"query": {"eo:cloud_cover": {"lt": max_cloud}},
}).encode()

req = urllib.request.Request(
"https://planetarycomputer.microsoft.com/api/stac/v1/search",
data=body, headers={"Content-Type": "application/json"},
)
features = json.load(urllib.request.urlopen(req, timeout=60))["features"]

# 반드시 직접 확인 — 검색 결과가 내 지점을 덮는다는 보장은 없다
covering = [f for f in features
if f["bbox"][0] <= lon <= f["bbox"][2]
and f["bbox"][1] <= lat <= f["bbox"][3]]
if not covering:
raise RuntimeError(f"{date_str} 에 해당 위치를 덮는 장면이 없습니다")

covering.sort(key=lambda f: f["properties"]["eo:cloud_cover"])
return covering[0]

이 함수를 실행하면 터미널에 이런 값이 찍힙니다.

2024: 2024-05-18 구름 0.00% 타일 52SDH 처리기준 05.10
2026: 2026-05-16 구름 8.71% 타일 52SDH 처리기준 05.12

API가 뭐라고 하는지 브라우저로 직접 보기​

프로그램: 브라우저 (Chrome)

검색 API는 URL 만 가지고 주소창에서도 열 수 있습니다. 아래 주소를 그대로 붙여보세요. 위 코드가 받은 것과 같은 JSON이 나옵니다.

https://planetarycomputer.microsoft.com/api/stac/v1/search?collections=sentinel-2-l2a&bbox=128.45,37.90,128.65,38.05&datetime=2024-05-18T00:00:00Z/2024-05-18T23:59:59Z&limit=2

처음 보면 무슨 소리인지 모릅니다. 그중 이 글에서 쓰는 다섯 개만 알면 됩니다.

키뜻이번 값
datetime촬영 시각2024-05-18
eo:cloud_cover타일 전체 구름 비율0.003637
s2:mgrs_tileMGRS 타일 번호52SDH
s2:processing_baseline방사 보정 버전05.10
bbox이 장면이 덮는 좌표 범위128.24129.11, 37.8638.85

그리고 장면 하나를 눈으로 보고 싶다면 아래 주소를 열면 그 위성사진이 지형 위에 펼쳐집니다.

https://planetarycomputer.microsoft.com/api/data/v1/item/map?collection=sentinel-2-l2a&item=S2B_MSIL2A_20240518T020649_R103_T52SDH_20240518T053353

브라우저에서 연 Planetary Computer의 장면 미리보기. 짙은 사각형이 52SDH 타일이고 설악산 해안선이 그 왼쪽 끝에 걸칩니다. 이 타일의 대부분은 바다이고 우리가 쓰는 3.8km는 그 안의 아주 작은 부분입니다.

여기서 중요한 게 하나 보입니다. 이 장면의 대부분은 바다입니다. 우리가 쓰는 대관령은 타일 끝에 딱 걸친 곳입니다. 그래서 앞서 말한 "검색 결과가 내 지점을 덮지 않는다"는 문제가 실제로 일어나는 겁니다.

eo:cloud_cover도 주의해서 읽으세요. 이건 타일 전체의 구름 비율이라, 내가 보는 3.8km 안에 구름이 없는지까지 알려주지 않습니다. 그건 SCL 밴드로 따로 확인해야 합니다.

Step 4. 밴드 하나를 내려받습니다​

프로그램: Python (urllib로 다운로드, Pillow로 열기) → 미리보기로 확인

밴드란 위성사진을 색상별로 분리한 층입니다. 필요한 밴드는 이름으로 불러옵니다.

밴드파장 영역이 글에서의 역할
B02청색진짜 색의 B 채널
B03녹색진짜 색의 G 채널
B04적색진짜 색의 R 채널, NDVI 계산값
B08근적외선사람 눈에는 안 보임, 식생에 반응이 큼

가장 중요한 부분: 아래의 min_value와 max_value를 반드시 고정해야 합니다.

REFLECT_MAX = 3000 # L2A 밴드의 반사율 표시 상한

def fetch_band(item_id, band, x, y, z, out_path):
url = (
f"https://planetarycomputer.microsoft.com/api/data/v1/item/tiles/WebMercatorQuad/{z}/{x}/{y}@1x.png"
f"?collection=sentinel-2-l2a&item={item_id}&assets={band}"
f"&min_value=0&max_value={REFLECT_MAX}"
)
req = urllib.request.Request(url, headers={"User-Agent": "tutorial"})
with urllib.request.urlopen(req, timeout=90) as resp:
data = resp.read()
with open(out_path, "wb") as f:
f.write(data)
return np.array(Image.open(out_path).convert("RGB")).astype(np.float32)[..., 0]

이 두 값을 빼면 밴드마다 밝기 늘림 범위가 제각각이 되어, 진짜 색이 검게 깨지거나 밴드 간 밝기 균형이 완전히 틀어집니다. 이 글 쓰는 동안 실제로 이 오류를 두 번 겪었습니다. 첫 번째는 인자 순서를 틀렸고, 두 번째는 이 값을 빠뜨린 경우였습니다.

내려받은 밴드를 미리보기로 열어봅니다​

프로그램: 미리보기 (Preview)

2024_B04.png를 열어보면 이런 화면이 나옵니다.

적색 밴드만 단독으로 연 화면. 밝기 보정 전이라 거의 검정으로 나옵니다. 이상한 오류가 아닙니다.

거의 검정인 게 정상입니다. 위성 밴드는 화면용 RGB가 아니라 반사율 숫자이고, 03000 범위를 0255로 압축한 결과 어두운 값이 많이 나오기 때문입니다. 이걸 그대로 화면으로 보여주면 안 됩니다. 그래서 색 조합이 필요합니다.

Step 5. 사람이 보는 색을 만듭니다​

프로그램: Python + matplotlib

밴드를 나란히 놓으면 색깔이 없습니다. 따라서 RGB 순서로 조합합니다. 이 순서가 위성 밴드 표기법의 표준입니다.

def truecolor(red, green, blue, lo=2, hi=98):
stack = np.dstack([red, green, blue])
low, high = np.percentile(stack, [lo, hi])
return np.clip((stack - low) / max(high - low, 1e-6) * 255, 0, 255).astype(np.uint8)

같은 밴드를 단독으로 연 화면과, RGB로 조합하고 밝기를 보정한 화면의 비교입니다.

lo, hi는 화면용 밝기 보정입니다. 계산에는 쓰지 않습니다. 화면용과 계산용을 분리하는 이 구분이 중요합니다. 밝기를 보정하면서 지수를 계산하면 숫자가 조작됩니다.

Step 6. 식생 지수(NDVI)를 계산합니다​

프로그램: Python + numpy로 계산, matplotlib으로 표시

NDVI는 "풀이 얼마나 살아 있는가"를 나타내는 값입니다. 식생은 근적외선을 강하게 반사하고 적색은 흡수하므로, 두 밴드의 차이가 커집니다.

def ndvi(nir, red):
return (nir - red) / np.maximum(nir + red, 1e-6)

B08과 B04로 계산한 NDVI 지도. 짙은 초록일수록 식생이 살아 있습니다.

값의 뜻은 대략 이렇습니다.

  • -1 ~ 0: 물, 그늘
  • 0 ~ 0.2: 노지, 맨땅
  • 0.2 ~ 0.4: 드문 풀
  • 0.4 ~ 0.7: 물이 적고 초록이 살아 있는 상태의 산림
  • 0.7 이상: 조밀한 수관

계산된 값은 색으로 바꿔서 봅니다.

ax.imshow(ndvi_value, cmap="YlGn", vmin=0, vmax=1)

cmap="YlGn"은 연노랑에서 짙은 초록으로 이어지는 색표입니다. vmin과 vmax를 0과 1로 고정해야 여러 시점·여러 지역 그림을 서로 비교할 수 있습니다. 비교의 전제가 시각화에서부터 무너지는 경우가 많습니다.

Step 7. 두 시점을 비교합니다 — 이게 제일 중요합니다​

프로그램: Python + matplotlib

여기까지는 욕심이 없습니다. 여기서 실수를하면 대번에 거짓말이 됩니다.

설악산 대관령을 2024년 5월과 2026년 5월 사진으로 겹쳤습니다. 2년 차이니 분명 변화가 보일 것 같고, 실제로 보입니다.

날짜 차이에 따른 같은 산의 NDVI 비교. 위는 날짜를 12일 어긋나게 고른 경우, 아래는 2일 차이로 맞췄을 때의 결과입니다. 같은 데이터와 같은 코드인데 결론이 완전히 다릅니다.

두 계산의 결과는 이렇습니다.

비교NDVI 평균차이 평균하락 0.15 초과
2024-05-18 vs 2026-05-06 (12일 차이)0.6465 → 0.5364-0.110125.34%
2024-05-18 vs 2026-05-16 (2일 차이)0.6465 → 0.6636+0.01710.09%

날짜를 12일 어긋나게 골랐을 때만, 타일의 4분의 1이 숲이 사라진 것처럼 나옵니다. 지도에 검은 선으로 표시된 면적이 25%인데, 실제 벌목이 아니었습니다.

그런데 2026년 장면 두 장의 구름 비율은 **1.02%와 8.71%**로 모두 구름이 거의 없는 장면입니다. 데이터가 나쁜 게 아닙니다. 5월 6일의 산과 5월 18일의 산은 다른 산이었던 겁니다. 한국에서 5월 중순은 새잎이 전개되는 시기이고, 5월 초과는 그보다 한창 이전입니다. 그 열흘짜리 차이를 빼면 식생 지수 전체가 기울어집니다.

그래서 2026년 5월 6일 장면이 더 좋게 보입니다. 구름이 적어서요. 날짜를 반대로 맞추면 오해를 그대로 갖고 갈 수 있는 구조입니다.

둘을 같은 계절로 맞춰 다시 계산하면 결론이 뒤집힙니다.

날짜를 2일 차이로 맞췄을 때의 결과. NDVI 차이 평균이 플러스 0.0171로 나오고, 0.15 이상 하락한 영역은 0.09%에 불과합니다.

  • 2년 평균이 오히려 0.017 올랐습니다. 이 3.8km 안에서는 뚜렷한 벌목 신호가 없습니다.
  • 검은 선으로 표시된 56개 화소(0.09%)는 흩어져 있고 덩어리를 이루지 않습니다. 벌목이라 부르기엔 근거가 없습니다.
  • "변화가 없다"가 하나의 결과입니다. 이 결론을 그대로 보고해야 합니다.

정리하면 비교할 때 지켜야 할 조건은 세 가지입니다.

  1. 연도 안에서 같은 계절로 맞춥니다. 연-월-일이 같은 시점이어야 합니다. 5월 초와 5월 중순은 다른 계절입니다.
  2. 구름은 SCL 밴드로 확인합니다. 타일 전체 구름 비율로 판단하면 안 됩니다.
  3. 처리 기준 버전을 기록합니다. 이번에는 2024년 장면이 05.10, 2026년 장면이 05.12였습니다. 같은 05.x 계열이라 영향이 작을 것으로 보이지만, 편향이 없다는 것은 이 글에서 검증하지 않았습니다. 장기 시계열을 만들 때는 반드시 같은 버전을 쓰거나, 버전을 맞춘 자료로 보정해야 합니다.

자주 막히는 곳​

화면에 나온 것실제로 잘못된 것해결
{"detail":"Tile(x=3161, y=7021, z=13) is outside bounds"}타일 번호를 z/y/x 순으로 넣었습니다순서를 z/x/y로 바꿉니다. lonlat_to_tile()이 x, y 순서로 돌려주므로 그대로 이어 붙이세요
{"detail":"assets must be defined either via expression or assets options."}assets를 쓰지 않았습니다&assets=B04 같은 항목을 최소 하나 넣습니다
{"detail":[{"type":"missing","loc":["query","collection"],"msg":"Field required"...}]}collection 또는 item이 없습니다URL에 collection=sentinel-2-l2a&item={장면ID}를 넣습니다
HTTP 500 응답rescale=minmax를 넣었습니다이 경로에서는 rescale가 500을 냅니다. min_value와 max_value를 쓰세요
이미지 파일이 1KB 남짓이고 거의 검정좌표가 타일 밖입니다, 또는 밝기 상한이 없습니다1. find_scene()의 bbox 검사를 통과했는지 확인 2. min_value/max_value가 있는지 확인
색이 진하게 초록으로 뜨거나 뭉개짐밝기 늘림 범위 밴드마다 다릅니다min_value=0, max_value=3000을 고정
구름이 있는데 왜 NDVI가 높게 나오죠타일 전체 구름 비율만 확인했습니다SCL 밴드로 내 3.8km 안의 구름을 확인
두 시점 비교에서 뜻밖의 큰 변화계절이 다릅니다연도 안에서 같은 계절로 날짜를 다시 고릅니다

전체 스크립트​

글에 나온 조각을 다 합치면 이 파일 하나로 끝납니다. 이 글의 모든 그림은 이 스크립트를 실행해 나온 결과입니다.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""위성이미지 완전 초보 튜토리얼 — 로그인 없이, 실제로 동작하는 최소 스크립트"""
import json
import math
import urllib.request

import numpy as np
from PIL import Image
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt

LON, LAT = 128.551, 37.952 # 설악산 대관령
ZOOM = 13 # 13이면 약 3.8km 정사각
REFLECT_MAX = 3000 # L2A 밴드 반사율 표시 상한


def lonlat_to_tile(lon, lat, z):
n = 2 ** z
x = int((lon + 180) / 360 * n)
r = math.radians(lat)
y = int((1 - math.log(math.tan(r) + 1 / math.cos(r)) / math.pi) / 2 * n)
return x, y


def find_scene(date_str, max_cloud=20):
day = f"{date_str}T00:00:00Z/{date_str}T23:59:59Z"
body = json.dumps({
"collections": ["sentinel-2-l2a"],
"bbox": [LON - 0.1, LAT - 0.05, LON + 0.1, LAT + 0.05],
"datetime": day,
"limit": 60,
"query": {"eo:cloud_cover": {"lt": max_cloud}},
}).encode()
req = urllib.request.Request(
"https://planetarycomputer.microsoft.com/api/stac/v1/search",
data=body, headers={"Content-Type": "application/json"},
)
features = json.load(urllib.request.urlopen(req, timeout=60))["features"]
# 검색 결과가 내 지점을 덮는다는 보장은 없다 — 반드시 직접 확인
covering = [f for f in features
if f["bbox"][0] <= LON <= f["bbox"][2] and f["bbox"][1] <= LAT <= f["bbox"][3]]
if not covering:
raise RuntimeError(f"{date_str} 에 해당 위치를 덮는 장면이 없습니다")
covering.sort(key=lambda f: f["properties"]["eo:cloud_cover"])
return covering[0]


def fetch_band(item_id, band, x, y, z, out_path):
url = (
f"https://planetarycomputer.microsoft.com/api/data/v1/item/tiles/WebMercatorQuad/{z}/{x}/{y}@1x.png"
f"?collection=sentinel-2-l2a&item={item_id}&assets={band}"
f"&min_value=0&max_value={REFLECT_MAX}"
)
req = urllib.request.Request(url, headers={"User-Agent": "tutorial"})
with urllib.request.urlopen(req, timeout=90) as resp:
data = resp.read()
with open(out_path, "wb") as f:
f.write(data)
return np.array(Image.open(out_path).convert("RGB")).astype(np.float32)[..., 0]


def ndvi(nir, red):
return (nir - red) / np.maximum(nir + red, 1e-6)


def truecolor(red, green, blue, lo=2, hi=98):
stack = np.dstack([red, green, blue])
low, high = np.percentile(stack, [lo, hi])
return np.clip((stack - low) / max(high - low, 1e-6) * 255, 0, 255).astype(np.uint8)


def main():
x, y = lonlat_to_tile(LON, LAT, ZOOM)
print(f"타일 번호: z{ZOOM} x{x} y{y}")

# 같은 계절로 맞춰 고른다. 날차를 크게 벌리면 계절 차이를 변화로 읽는다.
before = find_scene("2024-05-18")
after = find_scene("2026-05-16")

panels = {}
for label, scene in (("2024", before), ("2026", after)):
p = scene["properties"]
print(f"{label}: {p['datetime'][:10]} 구름 {p['eo:cloud_cover']:.2f}% "
f"타일 {p['s2:mgrs_tile']} 처리기준 {p['s2:processing_baseline']}")
item = scene["id"]
nir = fetch_band(item, "B08", x, y, ZOOM, f"{label}_B08.png")
red = fetch_band(item, "B04", x, y, ZOOM, f"{label}_B04.png")
grn = fetch_band(item, "B03", x, y, ZOOM, f"{label}_B03.png")
blu = fetch_band(item, "B02", x, y, ZOOM, f"{label}_B02.png")
panels[label] = {"ndvi": ndvi(nir, red), "rgb": truecolor(red, grn, blu)}

diff = panels["2026"]["ndvi"] - panels["2024"]["ndvi"]
print(f"\nNDVI 평균 2024: {panels['2024']['ndvi'].mean():.4f}")
print(f"NDVI 평균 2026: {panels['2026']['ndvi'].mean():.4f}")
print(f"차이 평균 : {diff.mean():+.4f}")
for t in (-0.05, -0.10, -0.15):
print(f" NDVI 하락 {abs(t):.2f} 초과 화소: {int((diff < t).sum())} ({(diff < t).mean() * 100:.2f}%)")

fig, ax = plt.subplots(1, 4, figsize=(22, 6))
ax[0].imshow(panels["2024"]["rgb"]); ax[0].set_title("2024-05-18 true color")
ax[1].imshow(panels["2026"]["rgb"]); ax[1].set_title("2026-05-16 true color")
im = ax[2].imshow(panels["2024"]["ndvi"], cmap="YlGn", vmin=0, vmax=1)
ax[2].set_title(f"NDVI 2024 (mean {panels['2024']['ndvi'].mean():.3f})")
im2 = ax[3].imshow(diff, cmap="RdYlGn", vmin=-0.35, vmax=0.35)
ax[3].contour(diff < -0.15, levels=[-0.149], colors="k", linewidths=1.0)
ax[3].set_title(f"NDVI change 2026-2024 (mean {diff.mean():+.3f})")
plt.colorbar(im, ax=ax[2], fraction=0.046)
plt.colorbar(im2, ax=ax[3], fraction=0.046)
for a in ax:
a.set_xticks([]); a.set_yticks([])
plt.suptitle("Seorak Daegwallyeong 37.952N 128.551E | Sentinel-2 L2A 10m | MGRS 52SDH | z13\n"
"Microsoft Planetary Computer, free, no login", fontsize=11)
plt.tight_layout()
plt.savefig("result.png", dpi=110)
print("\nresult.png 저장 완료")


if __name__ == "__main__":
main()

이 글에서 다루지 않은 것​

좋은 안내서는 어디까지가 유효 범위인지를 말해야 합니다. 다음은 직접 실행해 확인하지 않았으므로 확실한 방법이라고 적지 않습니다.

  • Google Earth Engine — 위성데이터 초보자에게 가장 많이 추천되지만 이 글의 경로와 완전히 다릅니다. 직접 해보고 별도로 정리했습니다.
  • Copernicus Data Space에서 원본 tif로 내려받기 — 해상도와 밴드가 더 좋지만 용량과 설정 난도가 큽니다. 이 글은 256×256 타일만 씁니다.
  • GDAL, rasterio — 큰 원본을 다루는 표준 도구입니다. 이 글에서 쓰지 않았습니다.
  • Landsat, Sentinel-1 레이더 데이터
  • 위성영상 분류 모델 학습 — 이 글은 지수 계산까지만 다룹니다.
  • 처리 기준 버전 차이(05.10 대 05.12)의 편향 정량화 — 버전이 다르고 그 사실만 확인했습니다. 편향이 없다는 근거는 없습니다.

포트폴리오·면접에서 이렇게 말해 보세요​

위성이미지를 다룬 프로젝트를 보여줄 때, 화면 수보다 비교 조건을 말하는 편이 훨씬 강합니다.

나쁜 예: "NDVI로 산림 변화를 추적했습니다."

좋은 예:

"같은 산을 2년 간격 위성사진으로 비교했습니다. 처음에 5월 초와 5월 중순을 비교해서 타일의 25%가 감소했다고 나왔는데, 알고 보니 열흘짜리 계절 차이였습니다. 같은 계절로 맞춰 다시 계산한 결과 0.09%였고, 이 지역에서는 뚜렷한 벌목 신호가 없다는 결론을 냈습니다. 다만 두 장면의 방사 보정 버전이 05.10과 05.12로 달라, 장기 시계열을 만들 때는 버전을 맞춰야 한다는 한계를 함께 기록했습니다."

이 답변에는 문제 발견, 검증, 오류 수정, 한계 인식이 다 들어 있습니다. 데이터 분석 면접에서 이 빈도를 말하는 사람은 많지 않습니다.

다음 액션​

  1. 위 스크립트를 그대로 한 번 돌려 result.png가 생기는지 확인합니다.
  2. LON, LAT만 자기 관심 지역으로 바꿔 다시 돌려봅니다.
  3. Step 7을 일부러 날짜를 12일 어긋나게 바꿔 돌리고, 25%대 수치가 나오는 것을 직접 확인합니다. 오해를 재현해 보는 경험이 설명보다 오래 남습니다.

다음 단계로는 지수만으로는 부족합니다. NDVI는 식생 유무를 말해주지만, 벌목인지 가뭄인지 병해인지 구분하지 못합니다. 거기서부터는 다중 시점 분석이나 분류 모델이 필요합니다.