기상레이더란 전파를 대기 중에 발사하여 강수입자에 부딪혀 산란되어 되돌아오는 신호를 수신·분석하여 강수구름의 위치, 강우강도, 이동속도 등을 탐지하는 기상관측장비입니다. 기상레이더는 1922년 테일러(Taylor)와 영(Young)이 지나가던 선박에 의해 신호가 변화하는 현상을 발견하면서 시작되어 제2차 세계대전을 통해 급속하게 발전하게 되었습니다. 적의 항공기나 선박을 탐지할 목적으로 개발된 군용 레이더는 제2차 세계대전이 끝난 후 일부가 기상용 레이더로 발전하였습니다.

기상레이더는 비(非)도플러 레이더(1세대), 단일편파 도플러 레이더(2세대), 이중편파 도플러 레이더(3세대) 순으로 발전해 왔습니다. 1세대 비도플러 레이더는 관측 지점에서 강수구름까지의 거리, 강수구름의 분포, 강수구름의 레이더반사도를 이용해 강수량을 추정할 수 있었습니다. 이후 1990년대에는 바람의 관측도 가능한 2세대 단일편파 도플러 레이더로 교체되었습니다.

레이더가 처음 발사한 전파와 강수입자에 반사되어 수신된 전파의 위상차를 이용하여 강수입자의 이동방향과 속도를 결정해 바람을 관측할 수 있게 되었습니다. 3세대 기상레이더는 이중편파 도플러 레이더로 수평으로 진동하는 전파와 수직으로 진동하는 전파를 동시에 발사하여 강수량을 추정할 뿐만 아니라, 비, 눈, 우박 등 강수형태를 구별할 수도 있습니다.

기상레이더는 대기중으로 전자기파를 대기중으로 발사하여 대기 중 강수입자로부터 반사되어 오는 신호를 수신하여 이로부터 강수입자의 강도와 위치를 측정하는 원격관측 장비 입나다.

- Radio : 고주파수의 전자기파(RF: Radio Frequency), 일반적으로 마이크로파를 지칭 - Detection : 물체의 존재를 탐지 - Ranging : 관측지점으로부터 목표물까지의 거리를 파악

우리나라의 기상레이더 관측은 1969년 11월 기상청이 국내 최초의 기상레이더를 서울 관악산에 설치하면서 시작되었습니다. 이후 백령도에서 제주도에 이르기까 지 전국에 9대의 현업용 기상레이더를 추가로 도입하여 총 10대로 구성된 기상레 이더 관측망을 운영하고 있습니다. 여기에 공항용 1대, 시험용 1대(레이더테스트 베드), 연구용 3대를 운영하고 있습니다.

우리나라와 주변에서 발달하는 집중호우, 태풍, 대설 등 위험기상 현상을 보다 신 속하고 정확하게 탐지할 수 있도록, 2014년 백령도를 시작으로 2019년까지첨단 성능을 가진 이중편파 기상레이더 관측망을 구축하여 운영 하고 있습니다.

기상레이더 정보는 다음 4단계를 거쳐 사용자에게 전달됩니다.

1. 관측 : 기상레이더가 대기 중으로 내보낸 전파는 비구름에 반사되어 일부가 약한 신호의 형태로 되돌아옵니다.
2. 신호처리 : 기상레이더는 수신된 신호를 신호처리기를 통해 변환하여 관측자료를 생산합니다.
3. 품질관리 : 관측자료로부터 기상현상이 아닌 정보를 걸러냅니다.
4. 정보 생산·제공: 품질관리된 자료를 활용하여 다양한 기상정보를 생산하여 제공합니다.

PPI자료는 레이더 볼륨자료 중 주어진 하나의 고도각에 대해서 관측된 자료 값을 말합니다. 기상청에서는 그 중 가장 낮은 고도각(PPI0)자료에 대해서 산출하고 있습니다. 레이더 관측 특성상 레이더로부터 멀어질수록 관측 고도가 높아지는 특징을 가지고 있습니다.

여러층의 PPI 관측자료 중 동일고도면(기상청에서는 보통 1.5km)의 자료를 이용하여 산출한 레이더 추정 강수량입니다. 동일한 고도의 강수량 분석에 유리합니다.

CMAX는 전체 고도각 자료 중 가장 강한 에코 부분의 자료를 산출한 자료입니다.

HSR은 레이더 불륨자료에서 지형에코, 차폐, 비(非)기상에코의 영향이 없으며 지면에 가장 가까운 고도각의 자료를 이용하여 산출한 레이더 추정 강수량입니다.

 사이트 자료를 관측 중첩지역에 대해 최댓값, 평균값, 근거리, 거리가중, 고도가중 등의 합성 기법에 따라 합성자료를 산출합니다. 

레이더 UF 자료는 연속된 Ray(Header + Data)를 반복하며 순차적으로 기록되어 있습니다.

1) Header 구조

① Mandatory Header : 위·경도 해발고도, 레이더명, 관측시각, sweep 번호, Ray 번호, 고도각, 방위각 등

② Optional Header, Local use header: 사용하지 않음 

③ Data Header : 자료종류와 각 자료의 입력 위치

④ Field Header : 빔폭, Nyquist velocity, 편파기능 유무,  Bin수, Scale factor 등

2) Data 구조

① 2 byte integer로 구성 

② 실제값 x Scale Factor (=100)

③ Null 값 = -32768

UF 레이더자료를 처리하기 위해서는 NASS TRMM Office에서 제공하는 RSL(Radar Software Library)가 필요하며 해당 라이브러리는 C언어용 라이브러리로 NASS 홈페이지에서 다운로드가 가능합니다. - 다운로드 페이지 : https://trmm-fc.gsfc.nasa.gov/trmm_gv/software/rsl/

다음은 RSL 라이브러리를 이용하여 레이더자료를 읽고 그림 그리는 예시입니다.

#include “rsl.h”

void main(int argc, char **argv)
{
Radar *radar;
radar = RSL_uf_to_radar(“radar.uf”);  // 레이더자료를 읽어 레이더 구조체로 저장 
RSL_load_refl_color_table();  //컬러테이블 정보를 읽음 
RSL_volume_to_gif(radar->v[CZ_INDEX], “cz_sweep”, 400, 400, 200.0)  //이미지 만들기

자세한 사용 매뉴얼은 홈페이지 또는 라이브러리 설치파일 “doc” 폴더에 있으니 참고하시면 됩니다.

NetCDF(Network Common Data Form) 포맷은 미국 기상연구 대학 연합(UCAR)에서 제공하는 포맷으로 배열지향적인 데이터를 쉽게 접근하고 활용할 수 있도록 개발한 포맷 형식입니다. 관련 라이브러리 및 매뉴얼은 https://unidata.ucar.ecu 에서 내려받아 활용할 수 있습니다. 홈페이지에서는 C, C++, 포트란, Java 라이브러리를 제공하고 있으며, 자제적으로 R, IDL, MATLAB, Ruby, Perl 등에서 데이터를 활용할 수 있는 인터페이스를 제공하고 있습니다.

레이더 사이트정보 NetCDF 구조는 다음과 같습니다.

netcdf RDR_BRI_RAW_202305241740 {
dimensions:
	time = 3240 ;
	n_points = 2534400 ;
	range = 960 ;
	sweep = 9 ;
	frequency = 1 ;
	string_length_sm = 10 ;
	string_length_md = 25 ;
variables:
	int volume_number ;
		volume_number:standard_name = "data_volume_index_number" ;
		volume_number:units = "1" ;
	double latitude ;
		latitude:standard_name = "latitude" ;
		latitude:units = "degrees_north" ;
	double longitude ;
		longitude:standard_name = "longitude" ;
		longitude:units = "degrees_east" ;
	double altitude ;
		altitude:standard_name = "altitude" ;
		altitude:units = "meters" ;
  ~~~ 중략

// global attributes:
		:Conventions = "Cf/Radial" ;
		:version = "1.3" ;
		:title = "vol2netcdf" ;
		:institution = "KMA" ;
		:references = "WRC" ;
		:source = "main.cpp" ;
		:history = "Phase offset" ;
		:comment = "2020.01.07" ;
		:instrument_name = "BRI" ;
		:site_name = "BRI" ;
		:scan_name = "BRI_240km_p2_vol" ;
		:scan_id = 1 ;
		:platform_is_mobile = "true" ;
		:ray_times_increase = "true" ;
		:field_names = "UH, UV, DBZH, DBZV, VELH, VELV, SWH, SWV, ZDR, PHIDP, RHOHV, NCPH, NCPV, SNRHC, SNRVC, CCOTH, CCORV" ;
		:time_coverage_start = "2023-05-24T08:40:01Z" ;
		:time_coverage_end = "2023-05-24T08:43:59Z" ;
}

기상청은 기상레이더를 이중편파레이더로 교체하면서 관측변수 추가(8종 → 16종) 및 정밀도 개선(8bit → 16bit)로 개선하기 위해 사이트자료 저장형식을 2019년부터 UF 포맷에서 NetCDF4으로 변경하여 운영하고 있습니다. UF 자료와 비교하여 저장되는 변수는 다음과 같습니다.

※ UF자료 중 PH, RH, DR, KD, HC는 이중편파레이더에만 저장되고 있음

레이더 합성자료는 기상청에서 정의한 이진파일 형태로 생산‧제공하고 있습니다. 합성자료 종류로는 HSR, CAPPI, PPI0, CMAX자료가 제공되고 있으며, 5분간격, 500m 수평해상도 자료로 제공되고 있습니다. 제공되는 레이더합성자료 자료구조는 다음과 같습니다.

struct RDR_CMP_HEAD {
  char   version;       // 포맷 버전
  short  ptype;         // 생산Product
  struct TIME_SS tm;    // 레이더 관측시각
  struct TIME_SS tm_in;// 합성자료 생성시각
  char   num_stn;       // 합성에 사용된 레이더 지점수
  char   map_code;      // 지도정보 코드
  char   map_etc;       // 기타 지도 정보코드(예비)
  short  nx;            // X축 격자점수 (1개 격자당 2 bytes로 저장함)
  short  ny;            // Y축 격자점수
  short  nz;            // Z축 격자점수
  short  dxy;           // 격자점간의 수평거리(m)
  short  dz;            // 격자점간의 수직거리(m) (nz=1이면 dz=0)
  short  z_min;         // nz > 1인 경우, 최저고도값(m) (nz <= 1이면 0)
  char   num_data;      // (nx*ny*nz)를 1개 자료블럭으로 했을때, 저장된 자료블럭수
  char   data_code[16];// 저장된 자료블력별 특성 코드
  char   etc[15];};
 };

struct RDR_CMP_STN_LIST {
  char  stn_cd[6];      // 레이더지점코드 
  struct TIME_SS tm;    // 레이더 관측시각
  struct TIME_SS tm_in; // 합성자료 생성시각
 };

struct TIME_SS {
  short  YY;    //년
  char   MM;   //월
  char   DD;    //일
  char   HH;    //시
  char   MI;    //분
  char   SS;    //초
};

레이더 합성자료 값이 반사도로 제공되는 경우 강우강도로 변환하여 활용해야 하는 경우가 있습니다. 기상청에서는 일반적으로 Z-R 관계식(Marshall-Palmer)을 활용하여 반사도를 강우강도로 변환합니다. 변환식은 다음과 같습니다(Z:반사도, R:,강우강도)

Z = 200 R1.6

반사도에서 강수강도로 상호변환하는 함수 예제코드(C언어)는 다음과 같습니다.

ZRa = 200;
 ZRb = 6;
// mode : 0(dBZ->강수량), 1(강수량->dBZ)
int dbz_rain_conv(float *dbz, float *rain, int mode)
{
  static int first = 0;
  static float za, zb;

  if (first == 0) {
    za = 0.1/var.ZRb;
    zb = log10(var.ZRa)/var.ZRb;
    first = 1;
  }
  if (mode == 0) {   
    *rain = *dbz*za - zb;
    *rain = pow(10.0, *rain);
  }
  else if (mode == 1) {
    *dbz = 10.0 * log10( var.ZRa * pow(*rain, var.ZRb) );
  }
  return 0;
}

• 사이트 : 기상자료개방포털 (data.kma.go.kr)​
• 메뉴 : 데이터 → 레이더 → 사이트, 합성

• 조건 : 자료형태, 기간, 지점, 시간

• 사이트 : 기상청 API허브 (apihub.kma.go.kr)
• 메뉴 : 레이더 → 레이더 강수량(HSR), 레이더 강수량, 레이더 원시자료 등

• 절차: 회원가입 및 로그인 → API 메뉴 이동 → API 활용신청 → API호출
• 호출 방법: 기상청_api허브_사용_안내서.pdf

소스코드: rdr_rain_notebook.ipynb

# API 요청을 위한 requests 라이브러리와 
# 분포도 표출을 위한 matplotlib 라이브러리,
# 이미지 및 점의 좌표계 변환을 위한 rasterio, pyproj 라이브러리와
# 지도를 표출하기 위한 folium 라이브러리를 설치합니다.
%pip install requests matplotlib rasterio pyproj folium

기상청 API허브에 회원가입 후 자신의 인증키를 사용하여 API 요청을 할 수 있습니다.

API 인증키를 코드상에 직접 남겨 사용하는 방식은 다른 사람에게 코드를 공유할때 인증키가 노출되는 문제가 생길 수 있습니다.

따라서 인증키를 환경변수로 등록하여 코드상에 노출되지 않도록 사용하는 것이 바람직합니다.

import os # 환경변수를 불러오는 표준 라이브러리

# 바람직하지 않은 방법
# my_api_key = 'MY_PRIVATE_API_AUTHENTICATE_KEY'

# 인증키를 'myapikey'라는 이름의 환경변수로 미리 등록해놓은 뒤 이를 불러와서 사용합니다.
my_api_key = os.environ['myapikey']

import requests # API 요청을 보내고 받는 라이브러리

# 레이더 합성자료 강수량 API url
api_url = 'https://apihub.kma.go.kr/api/typ04/url/rdr_cmp_file.php'

# 요청인자로 자료시간, 자료형태, 합성종류가 있습니다.

# 1. 자료시간을 년월일시분으로 입력합니다.
# 레이더 자료는 매 5분 단위로 입력할 수 있습니다.
time = '202208082000'

# 2. 자료 형태를 입력합니다.
# 이진파일(bin)과 이미지(img)를 입력할 수 있습니다.
data_type = 'bin'

# 3. 합성 종류를 입력합니다.
# 이진파일 요청의 경우 cpp, ppi. cmx, hsr, hsp를 입력할 수 있습니다.
cmp_type = 'hsp'

# API 요청인자들을 묶어 dictionary로 정의합니다.
api_parameters = {
    'tm': time,
    'data': data_type,
    'cmp': cmp_type,
    'authKey': my_api_key
}

# API허브에서는 이진 파일을 gz형식으로 압축하여 제공합니다.
file_name = f'RDR_CMP_{cmp_type.upper()}_PUB_{time}.bin.gz'

# 저장된 파일이 없는 경우에만 API를 요청해 파일을 다운로드합니다.
if not os.path.exists(file_name):
    # API 요청인자와 함께 API 요청
    response = requests.get(api_url, params=api_parameters)

    # 잘못된 응답을 받거나 짧은 에러메세지를 응답으로 받은 경우 에러 메세지 출력
    if response.status_code != 200 or len(response.content) < 500:
        raise requests.RequestException(
            response.content.decode('utf-8')
        )
    else:
        # 그 외의 올바른 응답에 대해서만 파일로 저장합니다.
        with open(file_name, 'wb') as f:
            f.write(response.content)
        print(f'{file_name} 파일 다운로드 완료')

다운로드 받은 레이더 합성자료에 관한 참고자료는 이곳에서 확인할 수 있습니다.

참고자료에 따르면 이진파일의 첫 1024byte는 헤더영역이며

나머지 byte는 2 byte(short int) * 2305(x좌표) * 2881(y좌표) 배열 영역으로 이루어져 있습니다.

배열 값 중 관측반경 바깥의 NULL값은 -30000 값을 사용하고 관측영역 내 비관측영역은 -25000 값을 사용하며

배열이 나타내는 실제 강수세기는 100을 나누어 사용해야 합니다.

자료의 설명에 맞게 gz파일의 압축을 풀어 강수세기 배열로 변환합니다.

import numpy as np  # 배열을 다루는 라이브러리
import gzip         # gz파일의 압축을 푸는 라이브러리

# 참고자료를 통한 배열의 x, y 크기를 정의합니다.
nx = 2305
ny = 2881

# 저장된 파일을 읽어 압축을 풉니다.
with open(file_name, 'rb') as f:
    decompressed_bytes = gzip.decompress(f.read())

# 파일의 1024번째 byte 이후로부터 short(np.int16) 자료형으로 읽어 배열을 정의합니다.
# 그리고 배열의 실제값은 정수값이 아닌 100을 나누어 계산한 소수값이므로
# 소수값을 담을 수 있도록 자료형을 np.int16(정수)에서 np.float32(부동소수점)으로 변환합니다.
# 마지막으로 배열의 형태를 2차원 x, y형태로 변환합니다.
# 이때, 좌표계의 x, y 순서와 배열의 행, 열은 순서가 반대인 것에 유의합니다.
rain_rate = np.frombuffer(
    decompressed_bytes,
    dtype=np.int16,
    offset=1024
).astype(np.float32).reshape(ny, nx)

# 배열의 관측반경 바깥의 영역을 나타내는 마스크를 정의합니다.
null_mask = (rain_rate <= -30000)

# 배열에서 관측반경 바깥의 영역값을 NaN값으로 변환합니다.
rain_rate[null_mask] = np.nan

# 마지막으로 배열의 값을 100으로 나누어 사용합니다.
rain_rate /= 100

# 배열의 행, 열 형태 확인
rain_rate.shape

강수 세기 분포도를 표출하기 위해 기상청 레이더 강수량 분포도에서 사용하는 색상표를 사용합니다.

# 분포도 표출을 위한 라이브러리
import matplotlib.pyplot as plt

# 분포도 표출을 위한 색상표
from matplotlib.colors import ListedColormap, BoundaryNorm, Normalize

# 레이더 강수량을 분포도로 표출하기 위해 25가지 색의 색상표를 RGB로 정의합니다.
colormap_rain = ListedColormap(np.array([
    [250, 250, 250], [0, 200, 255], [0, 155, 245], [0, 74, 245],                # 하늘색
    [0, 255, 0], [0, 190, 0], [0, 140, 0], [0, 90, 0],                          # 초록색
    [255, 255, 0], [255, 220, 31], [249, 205, 0], [224, 185, 0], [204, 170, 0], # 노랑색
    [255, 102, 0], [255, 50, 0], [210, 0, 0], [180, 0, 0],                      # 빨간색
    [224, 169, 255], [201, 105, 255], [179, 41, 255], [147, 0, 228],            # 보라색
    [179, 180, 222], [76, 78, 177], [0, 3, 144], [51, 51, 51]                   # 파란색
]) / 255)

# 색상표에서 NaN값은 투명한 색상(RGBA 중 A가 0)을 나타내도록 합니다.
colormap_rain.set_bad([0, 0, 0, 0])

# 색상표의 각각의 색이 나타내는 범위를 정한다.
bounds = np.array([
    0, 0.1, 0.5, 1, # 하늘색
    2, 3, 4, 5,     # 초록색
    6, 7, 8, 9, 10, # 노랑색
    15, 20, 25, 30, # 빨간색
    40, 50, 60, 70, # 보라색
    90, 110, 150    # 파란색
])

# 색상 개수와 값 범위를 서로 맞춥니다.
norm = BoundaryNorm(boundaries=bounds, ncolors=len(colormap_rain.colors))

# 색상표를 표출합니다.
colormap_rain

# 분포도에서 색상표 위치, 크기를 조절하는 함수
from mpl_toolkits.axes_grid1 import make_axes_locatable

# 강수세기 배열 각각의 값을 색상표의 색이 나타내는 인덱스로 매핑시킵니다.
# 매핑된 각 인덱스 범위(0 ~ 25)를 정규화(0 ~ 1)합니다.
# 정규화를 통해 잃게된 투명색의 범위를 다시 지정합니다.
# 마지막으로 색상표의 색상으로 배열을 매핑합니다.
# 이때 배열의 형태를 (행 개수, 열 개수) 에서 (행 개수, 열 개수, RGBA)로 변환합니다.
colored_array = BoundaryNorm(
    boundaries=bounds,
    ncolors=len(colormap_rain.colors)
)(rain_rate)
colored_array = Normalize(
    0, len(colormap_rain.colors)
)(colored_array)
colored_array[null_mask] = np.nan
colored_array = (colormap_rain(colored_array) * 255).astype(np.uint8)

# 분포도에 표시되는 범위 값을 2칸 간격으로 구합니다.
ticks = bounds[:]

# 분포도를 그릴 크기를 지정합니다.
fig, ax = plt.subplots(1, 1, figsize=(5, 5))

# 분포도의 제목을 지정합니다.
ax.set_title('Colored array')

# 배경색을 회색으로 지정합니다.
ax.set_facecolor('#cccccc')

# 색상표의 크기와 위치를 조절합니다.
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0)

# 배열을 해당 크기에 맞춰 그립니다.
im = ax.imshow(colored_array, origin='lower', cmap=colormap_rain, norm=norm)

# 색상표에 표시될 글자 크기 및 제목을 설정합니다.
cbar = fig.colorbar(im, cax=cax, ticks=ticks)
cbar.ax.tick_params(labelsize=8)
cbar.ax.set_title('mm/h', fontsize=8)

# 그린 이미지를 표출합니다.
plt.show()

위와 같은 이미지를 동적인 지도 위에 표출하기 전, 지도와 레이더 자료(이미지)의 좌표계(Coordinate Reference System, CRS)를 맞추어야 합니다.

현재 동적인 지도가 사용하는 좌표계는 웹 상에서 지도를 표출할 때 많이 사용되는 EPSG:3857입니다.

참고자료에 따르면 레이더 자료가 사용하는 좌표계는 Lambert Conformal Conic Projection(LCC)으로 한반도를 중심으로 표현하기 적합해 여러 데이터에서도 사용된다고 합니다.

레이더 자료를 EPSG:3857로 변환하여 지도 위에 표출합니다.

# 좌표계 변환을 위한 변환 행렬
from rasterio.transform import Affine

# 좌표계 변환을 위해 필요한 함수
from rasterio.warp import calculate_default_transform, reproject, Resampling

# 레이더 자료의 너비와 높이, 중심점 좌표와 공간 해상도를 m단위로 정의합니다.
source_width = nx
source_height = ny
source_center_x = 1121
source_center_y = 1681
source_resolution = 500

# 변환 전의 좌표계를 proj.4 형태의 문자열로 정의합니다.
# 투영법(Projection)은 LCC이며 LCC 좌표계 정의에 필요한 위도 선 2개(lat_1, lat_2)와
# 중심 위도, 경도(lat_0, lon_0)을 정의합니다.
# 좌표계의 x, y는 각각 오른쪽, 위로 증가하는 방향을 가지고(False easting, False northing)
# 지구 타원체를 WGS84로 정의합니다.
# 마지막으로 좌표계가 사용하는 단위는 미터(m)로 정의합니다.
source_crs = "+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=126 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"

# 이미지(배열)의 행과 열이 LCC 좌표계에서의 좌표로 변환되기 위한 변환 행렬(Affine Transform Matrix)을 정의합니다.
# 변환 행렬은 먼저 이미지의 중심점을 (0, 0) 위치로 이동시킨 뒤(Translation)
# 각 이미지 픽셀이 실제로 500m를 나타내므로 크기를 x, y로 500배 키워주는(Scale) 역할을 합니다.
# 이때, 행렬의 곱하는 순서에 유의합니다.(AB != BA)
source_transform = Affine.scale(source_resolution, source_resolution) * Affine.translation(-source_center_x, -source_center_y)

# 변환 행렬을 거친 이미지가 나타내는 경계를 정의합니다.
source_bounds = {
    'left': -source_center_x * source_resolution,
    'bottom': (source_height - source_center_y) * source_resolution,
    'right': (source_width - source_center_x) * source_resolution,
    'top': -source_center_y * source_resolution
}

# 변환 후 이미지의 변환 행렬과 너비와 높이를 계산합니다.
dest_transform, dest_width, dest_height = calculate_default_transform(
    src_crs=source_crs,
    dst_crs='EPSG:3857',
    width=source_width,
    height=source_height,
    **source_bounds,
)

# 변환 후의 이미지가 담길 비어있는 배열을 정의합니다.
converted_array = np.ones((dest_height, dest_width, 4), dtype=np.uint8)

# RGBA 각 채널에 대해 좌표계 변환을 수행합니다.
# 사용하는 resampling 기법으로 가까운 값을 선택하는 nearest를 선택합니다.
for i in range(4):
    reproject(
        source=colored_array[:, :, i],
        destination=converted_array[:, :, i],
        src_transform=source_transform,
        src_crs=source_crs,
        dst_transform=dest_transform,
        dst_crs='EPSG:3857',
        resampling=Resampling.nearest,
    )

# 변환된 이미지를 분포도로 표출합니다.
# 이전 분포도와 다르게 원점의 위치, 너비, 높이가 달라졌음을 확인합니다.
fig, ax = plt.subplots(1, 1, figsize=(5, 5))

# 분포도의 제목을 지정합니다.
ax.set_title('Converted array')

# 배경색을 회색으로 지정합니다.
ax.set_facecolor('#cccccc')

# 색상표의 크기와 위치를 조절합니다.
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0)

# 배열을 표출하며 정의한 색상표와 범위를 지정합니다.
im = ax.imshow(converted_array, cmap=colormap_rain, norm=norm)

# 색상표에 표시될 글자 크기 및 제목을 설정합니다.
cbar = fig.colorbar(im, cax=cax, ticks=ticks)
cbar.ax.tick_params(labelsize=8)
cbar.ax.set_title('mm/h', fontsize=8)

# 그린 이미지를 표출합니다.
plt.show()

마지막으로 변환된 이미지를 동적인 지도 위에 표출합니다.

import folium   # 지도로 표출하기 위한 라이브러리
import branca.colormap as cm    # 지도에 색상표를 표출하기 위한 색상표 모듈

# 지도에 마우스 위치를 표출하는 플러그인
from folium.plugins import MousePosition

# 점끼리의 좌표계 변환을 위한 transformer
from pyproj.transformer import Transformer

# folium을 이용한 좌표는 경도, 위도만 입력 가능하므로
# EPSG:3857 좌표계에서 EPSG:4326 좌표계로의 변환하는 transformer를 정의합니다.
degree_tranformer = Transformer.from_crs('EPSG:3857', 'EPSG:4326')

# 너비, 높이 크기가 800인 지도 생성
fig = folium.Figure(width=800, height=800)

# 지도의 중심과 나타낼 영역을 위·경도로 입력합니다.
map = folium.Map(
    location=[38, 126],
    zoom_start=6,
    min_zoom=6,
    min_lat=28,
    max_lat=43,
    min_lon=116,
    max_lon=135,
    max_bounds=True,
).add_to(fig)

# 지도에 색상표를 추가합니다.
map.add_child(cm.StepColormap(
    [tuple(i) for i in colormap_rain.colors],
    vmin=bounds[0], vmax=bounds[-1], tick_labels=[], caption='mm/h'
))

# 변환된 이미지를 지도 위에 겹쳐 그립니다.
# 불투명도를 0.4로 정하고
# 이미지가 그려질 범위를 경도, 위도로 입력합니다.
folium.raster_layers.ImageOverlay(
    image=converted_array,
    name='rain_rate',
    opacity=0.4,
    bounds=[
        degree_tranformer.transform(*dest_transform.__mul__((0, dest_height))),
        degree_tranformer.transform(*dest_transform.__mul__((dest_width, 0)))
    ],
    zindex=1,
).add_to(map)

# 지도에 마우스 위치를 표출하는 플러그인과 layer를 조절하는 버튼을 추가합니다.
MousePosition().add_to(map)
folium.LayerControl().add_to(map)

# 지도를 표출합니다.
fig

# 지도 상에서 원하는 점의 위·경도를 알아낸뒤, 그 값을 인덱싱하여 가져올 수 있습니다.
# 예를 들어, 부산 해운대의 위·경도가 38.19192, 128.59634 이라면
# 그 지점의 좌표를 다시 LCC 좌표계로 변환하여 그 값을 가져올 수 있습니다.

# 위·경도 좌표계에서 LCC 좌표계로 변환하는 transformer를 정의합니다.
lcc_transformer = Transformer.from_crs('EPSG:4326', source_crs)

# 원하는 지점의 위·경도를 입력합니다.
target_lat = 38.19192
target_lon = 128.59634

# 위·경도 값을 LCC 좌표계로 변환합니다.
# source_transform 연산의 반대이므로 역행렬을 곱합니다.
index_col, index_row = source_transform.__invert__().__mul__(
    lcc_transformer.transform(target_lat, target_lon)
)

# 지점의 강수세기를 출력합니다.
print(
    f"위도 {target_lat}, 경도 {target_lon} 지점의 {time} 강수 세기는 "
    f"{rain_rate[round(index_row), round(index_col)]:.02f}mm/h 입니다."
)