Try to draw the global rain fall based on JAXA GSMaP data～JAXA GSMaPを利用して全球の降雨量を描画してみます
This TimeもGeminiにAdviseをGiveしてもらいながら、DigitalWorldをExploreします。
今回は、GlobalなWaterResourceについてMeなりにTreasureHuntingしたいとSeemします。
For the first time in a while, I try to run a Python code. In long time, scince I work on my work of My Biz, I couldn't study Python code.
※GSMaP data by Japan Aerospace Exploration Agency (JAXA)"
※本論文にて使用したGSMaPデータは、宇宙航空研究開発機構より提供を受けました。
As reference, GSMaPはGlobal Satellite Mapping Precipitationの略だそうです。
GSMaPのHomePage(https://sharaku.eorc.jaxa.jp/GSMaP/index_j.htm)はLikeBelowのImageです。Upper Left SideにあるOrangeのButtonをPushするとUserRegisterPageにTranssionします。
UserRegisterPageはLikeBelowです。
Nuturallyですが、Required ItemsをInputして、Registerします。RegisterがCompleteすると、MaybeでGSMaPのFTP(File Transfer ProtcolというFilesを転送するためのProtcolSyatem)に関するFTP_HOSTのAdressやFTP_USERのID、FTP＿PASSがByE-MailでSendされてきたとRememberしています。
ちなみにFTP、とっても懐かしいです。わたしが大学生の時代からあり、LargeSizeFilesをTransferするときにUsedしていました。
ここからはGeminiとALotOfTimesもCatchBallしながら、GoAheadしていきました。
FTP_HOSTなどをGetすると、Geminiが早速Precipitation(=A Volume of RainFall)のDownloadCodeをOfferしてくれたのですが、そのCodeはFTPの中にあるFileのPathやNameをCorrectにRecognizeできておらず、Don'tWorkだったので、まずここでCatchBallがStartです。
ちなみにGeminiからはFTP利用にあたって汎用のFreeAppをUseすることをRecommendedされましたが、AppをInstallするのがHave a troubleだったので、PythonでCodeをWriteしてもらいました。
そこで、Geminiさんから「At FirstはFTP_ServerのContentsやDirectryStructureをReviewしてください」とのSugestionとCodeForItをOfferしてもらいました。
import ftplib
# --- License and Disclaimer ---
# Copyright (c) 2026 [Author Name]
# Released under the MIT License.
# 【免責事項】本スクリプトによるFTPアクセスおよび実行結果について、開発者は一切の責任を負いません。
# ------------------------------
#FTP_HOST = "XXXX.YYYY.jaxa.jp"
#FTP_USER = "YOUR_USERNAME" # JAXAから発行されたユーザーIDを入力
#FTP_PASS = "YOUR_PASSWORD" # JAXAから発行されたパスワードを入力
# メールに記載されていた確実な上位階層を指定
TARGET_DIR = "/standard/v8"
def explore_directory():
try:
with ftplib.FTP(FTP_HOST) as ftp:
ftp.login(FTP_USER, FTP_PASS)
ftp.cwd(TARGET_DIR)
print(f"=== Directory contents of {TARGET_DIR} ===")
# フォルダやファイルの一覧を出力
ftp.retrlines('LIST')
except ftplib.error_perm as e:
print(f"FTP Error: {e}")
except Exception as e:
print(f"Error: {e}")
if __name__ == "__main__":
explore_directory()
By the Way、MarkDownしている
#FTP_HOST = "XXXX.YYYY.jaxa.jp"
#FTP_USER = "YOUR_USERNAME" # JAXAから発行されたユーザーIDを入力
#FTP_PASS = "YOUR_PASSWORD" # JAXAから発行されたパスワードを入力
のPartは、SeparatelyにCodingしておくNecesarryがあります。
このCodeのOutputResultはBelowです。Briefingすると、JAXAのGSMaPのDATAを収めているFTP_Serverの\standard\v8というDirectoryのInsideのStructureをGoLookしているようです。
=== Directory contents of /standard/v8 ===
-rw-r--r-- 1 411 trmm 6373 Oct 9 2024 GSMaP_MVK_RNL_HISTORY.txt
-rw-r--r-- 1 411 trmm 6101 Oct 8 2024 README.first.txt
drwxr-xr-x 4 411 trmm 4096 Dec 8 2021 daily
drwxr-xr-x 4 411 trmm 4096 Dec 8 2021 daily_G
drwxr-xr-x 2 411 trmm 4096 Oct 9 2024 doc
drwxr-xr-x 31 411 trmm 4096 Jan 5 2026 hourly
drwxr-xr-x 31 411 trmm 4096 Jan 5 2026 hourly_G
drwxr-xr-x 31 411 trmm 4096 Feb 3 2026 monthly
drwxr-xr-x 31 411 trmm 4096 Feb 3 2026 monthly_G
drwxr-xr-x 31 411 trmm 4096 Jan 5 2026 netcdf
drwxr-xr-x 31 411 trmm 4096 Jan 5 2026 reliability
drwxr-xr-x 2 411 trmm 4096 Jul 10 2023 sample
drwxr-xr-x 31 411 trmm 4096 Jan 5 2026 sateinfo
drwxr-xr-x 31 411 trmm 4096 Jan 5 2026 timeinfo
drwxr-xr-x 5 411 trmm 4096 Dec 8 2021 txt
ThisOutputをGeminiにGiveすると、BelowのようなAnswerをいただきました。
- 仕様確認（README）: README.first.txt にデータ形式（バイナリかNetCDFか）、欠損値の値、解像度などの重要仕様が記載されています。
- ターゲットディレクトリ: 地上の雨量計データで誤差補正された高精度な月間データが格納されている monthly_G （Gauge-adjusted）が最適です。また、当初想定していたNetCDF形式のデータが netcdf ディレクトリ内にまとめられている可能性もあります。
- 次は
README.first.txtの中身をテキストとして読み込みつつ、最も有望なmonthly_Gフォルダのさらに下の階層（年号フォルダが存在するかどうか）を確認します。以下のコードを実行してみてください。
Jupyter Notebook等の前段セルで認証情報が定義されている前提に基づき、ヘッダー情報を統合し、認証部分をコメントアウトした探索コードです。このセルを実行することで、
README.first.txtの全文とmonthly_G直下の年号ディレクトリ一覧が出力されます。import ftplib
# ==============================================================================
# [Feature]
# JAXA GSMaPのFTPサーバーにアクセスし、データ仕様書(README.first.txt)のテキスト内容を
# コンソールに直接出力するとともに、雨量計補正版の月間データが格納されている
# 「monthly_G」ディレクトリの内部構造(年号フォルダ等)を探索・確認します。
#
# [License & Disclaimer]
# Copyright (c) 2026 [Author Name] / Released under the MIT License.
# 【免責事項】本スクリプトによるFTPアクセスおよび実行結果について、開発者は一切の責任を負いません。
# ==============================================================================
#FTP_HOST = "hokusai.eorc.jaxa.jp"
#FTP_USER = "YOUR_USERNAME" # JAXAから発行されたユーザーID
#FTP_PASS = "YOUR_PASSWORD" # JAXAから発行されたパスワード
def read_readme_and_explore_monthly():
with ftplib.FTP(FTP_HOST) as ftp:
ftp.login(FTP_USER, FTP_PASS)
print("=== Contents of README.first.txt ===")
ftp.cwd("/standard/v8")
try:
# READMEの中身を直接コンソールに出力
ftp.retrlines('RETR README.first.txt')
except Exception as e:
print(f"Error reading README: {e}")
print("\n=== Directory contents of /standard/v8/monthly_G ===")
try:
# 月間データ(雨量計補正版)のディレクトリの中身を確認
ftp.cwd("/standard/v8/monthly_G")
ftp.retrlines('LIST')
except Exception as e:
print(f"Error accessing monthly_G: {e}")
if __name__ == "__main__":
read_readme_and_explore_monthly()
Outputは Quite LongですBut、All OutputをCopy＆A¥Pasteです。
=== Contents of README.first.txt ===
README.first (GSMaP_MVK & GSMaP_Gauge Ver.8)
********************************************************************************
* Documents
********************************************************************************
Readme (this file): /standard/v8/README.first.txt
Version up history: /standard/v8/GSMaP_MVK_RNL_HISTORY.txt
Caveat for algorithm: /standard/v8/doc/CaveatForUseOfGPM-GSMaP_Ver.8.pdf
Data Format Document: /standard/v8/doc/DataFormatDescription_MVK_RNL_v8.0000.pdf
User's Guide: /standard/v8/doc/USERS_GUIDE_MVK_RNL_v8.0000.pdf
********************************************************************************
* Archive Period and Products
********************************************************************************
Since January 1, 1998
Global Satellite Mapping of Precipitation Microwave-IR Combined Product (GSMaP_MVK)
Gauge-calibrated Rainfall Product (GSMaP_Gauge)
********************************************************************************
* Data Directory (ver.8) and File Naming Rule
********************************************************************************
Hourly Rain Rate data;
Directory: /standard/v8/hourly/YYYY/MM/DD/
File Name: gsmap_mvk.YYYYMMDD.HHNN.vP.RSKI.J.dat.gz (GSMaP_MVK)
Hourly Gauge-calibrated Rain Rate data;
File Name: gsmap_gauge.YYYYMMDD.HHNN.vP.RSKI.J.dat.gz (GSMaP_Gauge)
Satellite Information Flag;
Directory: /standard/v8/sateinfo/YYYY/MM/DD/
File Name: gsmap_mvk.YYYYMMDD.HHNN.vP.RSKI.J.sateinfo.dat.gz (GSMaP_MVK)
Observation Time Flag;
Directory: /standard/v8/timeinfo/YYYY/MM/DD/
File Name: gsmap_mvk.YYYYMMDD.HHNN.vP.RSKI.J.timeinfo.dat.gz (GSMaP_MVK)
Reliability Flag;
Directory: /standard/v8/reliability/YYYY/MM/DD/
File Name: gsmap_mvk.YYYYMMDD.HHNN.vP.RSKI.J.reliability.dat.gz (GSMaP_MVK)
Hourly Rain Rate & major flag in NetCDF;
Directory: /standard/v8/netcdf/YYYY/MM/DD/
File Name: gsmap_mvk.YYYYMMDD.HHNN.vP.RSKI.J.nc
Hourly Rain Rate & major flag in HDF;
Directory: /HDF/standard/v8/Hourly/YYYY/MM/DD/
File Name: GPMMRG_MAP_YYMMDDHHNN_H_L3S_MCH_VVV.h5
Daily Averaged Rain Data (00Z-23Z averaged);
Directory: /standard/v8/daily/00Z-23Z/YYYYMM/
File Name: gsmap_mvk.YYYYMMDD.0.1d.daily.00Z-23Z.vP.RSKI.J.dat.gz (GSMaP_MVK)
Daily Averaged Rain Data (p12Z-11Z averaged);
Directory: /standard/v8/daily/p12Z-11Z/YYYYMM/
File Name: gsmap_mvk.YYYYMMDD.0.1d.daily.p12Z-11Z.vP.RSKI.J.dat.gz (GSMaP_MVK)
Daily Averaged Gauge-calibrated Rain Data (00Z-23Z averaged);
Directory: /standard/v8/daily_G/00Z-23Z/YYYYMM/
File Name: gsmap_gauge.YYYYMMDD.0.1d.daily.00Z-23Z.vP.RSKI.J.dat.gz (GSMaP_Gauge)
Daily Averaged Gauge-calibrated Rain Data (p12Z-11Z averaged);
Directory: /standard/v8/daily_G/p12Z-11Z/YYYYMM/
File Name: gsmap_gauge.YYYYMMDD.0.1d.daily.p12Z-11Z.vP.RSKI.J.dat.gz (GSMaP_Gauge)
Hourly Rain & Gauge-calibrated Rain Data in Text Format;
Directory: /standard/v8/txt/hourly/XX_ZZZZZZ/YYYY/MM/DD/
File Name: gsmap_mvk_vPRSKIJ_YYYYMMDD_HH00_XX_ZZZZZZ.csv.zip (GSMaP_MVK & GSMaP_Gauge)
Daily Averaged Rain & Gauge-calibrated Rain Data in Text Format (00Z-23Z averaged);
Directory: /standard/v8/txt/daily/00Z-23Z/XX_ZZZZZZ/YYYY/MM/
File Name: gsmap_mvk_vPRSKIJ_YYYYMMDD_daily_00Z-23Z_XX_ZZZZZZ.csv.zip (GSMaP_MVK & GSMaP_Gauge)
Daily Averaged Rain & Gauge-calibrated Rain Data in Text Format (p12Z-11Z averaged);
Directory: /standard/v8/txt/daily/p12Z-11Z/XX_ZZZZZZ/YYYY/MM/
File Name: gsmap_mvk_vPRSKIJ_YYYYMMDD_daily_p12Z-11Z_XX_ZZZZZZ.csv.zip (GSMaP_MVK & GSMaP_Gauge)
Daily Averaged Rain & major flag in HDF;
Directory: /HDF/standard/v8/Daily/YYYY/MM
File Name: GPMMRG_MAP_YYMMDD_D_L3S_MCD_VVV.h5
Monthly Averaged Rain Data;
Directory: /standard/v8/monthly/YYYY/
File Name: gsmap_mvk.YYYYMM.0.1d.monthly.vP.RSKI.J.dat.gz (GSMaP_MVK)
Monthly Averaged Gauge-calibrated Rain Data;
Directory: /standard/v8/monthly_G/YYYY/
File Name: gsmap_gauge.YYYYMM.0.1d.monthly.vP.RSKI.J.dat.gz (GSMaP_Gauge)
Monthly Averaged Rain & Gauge-calibrated Rain Data in Text Format;
Directory: /standard/v8/txt/monthly/XX_ZZZZZZ/YYYY/
File Name: gsmap_mvk_vPRSKIJ_YYYYMM_monthly_XX_ZZZZZZ.csv.zip (GSMaP_MVK & GSMaP_Gauge)
Monthly Averaged Rain & major flag in HDF;
Directory: /HDF/standard/v8/Monthly/YYYY
File Name: GPMMRG_MAP_YYMM_M_L3S_MCM_VVV.h5
where,
YYYY: 4-digit year
YY: 2-digit year
MM: 2-digit month
DD: 2-digit day
HH: 2-digit hour
NN: 2-digit minute (currently fixed as 00)
P: Algorithm version
R: Version of microwave imager algorithm (reset when P is updated)
S: Version of microwave sounder algorithm (reset when P is updated)
K: Version of microwave imager/sounder algorithm (reset when P is updated)
I: Version of microwave-IR combined algorithm
J: Inclement number of reprocessing
XX_ZZZZZZ: 9-digit area name
VVV: 3-digit product version.
********************************************************************************
* Sample code
********************************************************************************
Directory: /standard/v8/sample/
Fortran Sample Code: read_GSMaP_MVK_0.1deg.v8.f
IDL Sample Code: GSMaP_MVK_sample.v8.pro
Python sample code: readGSMaP_MVK_netcdf.py
GrADS Control File:
Hourly Rainfall Data GSMaP_MVK.hourly.rain.v8.ctl
Satellite Information Flag GSMaP_MVK.hourly.sat.v8.ctl
Observation Time Flag GSMaP_MVK.hourly.time.v8.ctl
Reliability Flag GSMaP_MVK.hourly.reliability.v8.ctl
Daily (00Z-23Z avaraged) GSMaP_MVK.daily.00Z-23Z.v8.ctl
Daily (p12Z-11Z averaged) GSMaP_MVK.daily.p12Z-11Z.v8.ctl
Monthly GSMaP_MVK.monthly.v8.ctl
********************************************************************************
=== Directory contents of /standard/v8/monthly_G ===
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 1998
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 1999
drwxr-xr-x 2 411 trmm 4096 Sep 14 2023 2000
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2001
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2002
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2003
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2004
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2005
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2006
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2007
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2008
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2009
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2010
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2011
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2012
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2013
drwxr-xr-x 2 411 trmm 4096 Jul 9 2023 2014
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 2015
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 2016
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 2017
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 2018
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 2019
drwxr-xr-x 2 411 trmm 4096 Sep 7 2023 2020
drwxr-xr-x 2 411 trmm 4096 Sep 14 2023 2021
drwxr-xr-x 2 411 trmm 4096 Mar 13 2023 2022
drwxr-xr-x 2 411 trmm 4096 Jan 3 2024 2023
drwxr-xr-x 2 411 trmm 4096 Jan 3 2025 2024
drwxr-xr-x 2 411 trmm 4096 Jan 3 2026 2025
drwxr-xr-x 2 411 trmm 4096 Sep 3 03:56 2026
ThisOutputをGeminiにGiveして、AdviseをOfferしてもらいました。
出力結果の共有ありがとうございます！ READMEの仕様と、
monthly_G 直下に 1998 から 2026 までの年号ディレクトリが揃っていることが完全に確認できました。READMEの記載から、今回のデータに関する重要な仕様が2つ判明しています。
ファイル形式: NetCDFではなく、GZIP圧縮されたRawバイナリ形式（
*.dat.gz）であること。データの単位: GSMaPの月間データは通常「mm/hr（1時間あたりの平均降水量）」で格納されているため、累計を出すには「データ × 24時間 × その月の日数」を掛けて月間総雨量（mm/month）に換算してから足し合わせる必要があること。
これらを踏まえ、FTPから dat.gz ファイルをダウンロードし、NumPyでバイナリを展開して年間累計降水量を計算し、GeoTIFFとして書き出す一連のコードを作成しました。
# ==============================================================================
# [Feature]
# JAXA GSMaP(monthly_G)のFTPサーバーから、過去10年間(2016〜2025年)の月間降水量データを
# 自動ダウンロードします。全120ヶ月分のデータ(mm/hr)を月ごとの日数を用いて総降水量に換算・累計し、
# 最後に10で割ることで「過去10年の平均年間降水量マップ」をGeoTIFFとして出力します。
#
# [License & Disclaimer]
# Copyright (c) 2026 [Author Name] / Released under the MIT License.
# 【免責事項】本スクリプトによるデータ取得および算出結果は衛星観測に基づく推計値であり、
# その正確性や特定目的への適合性を保証するものではありません。
# 開発者は本コードの利用により生じた一切の損害について責任を負いません。
# ==============================================================================
import ftplib
import os
import gzip
import calendar
import numpy as np
import rasterio
from rasterio.transform import from_origin
# FTP_HOST = "hokusai.eorc.jaxa.jp"
# FTP_USER = "YOUR_USERNAME"
# FTP_PASS = "YOUR_PASSWORD"
START_YEAR = 2016
END_YEAR = 2025 # 2016〜2025の10年間
BASE_DIR = "./gsmap_monthly_G_10years"
OUTPUT_TIFF = "gsmap_10yr_avg_annual_precipitation.tif"
def download_10years_data():
os.makedirs(BASE_DIR, exist_ok=True)
with ftplib.FTP(FTP_HOST) as ftp:
ftp.login(FTP_USER, FTP_PASS)
for year in range(START_YEAR, END_YEAR + 1):
ftp_dir = f"/standard/v8/monthly_G/{year}/"
try:
ftp.cwd(ftp_dir)
files = ftp.nlst("*.dat.gz")
for file in files:
local_path = os.path.join(BASE_DIR, file)
if not os.path.exists(local_path):
print(f"Downloading {file}...")
with open(local_path, "wb") as f:
ftp.retrbinary(f"RETR {file}", f.write)
except Exception as e:
print(f"Error accessing directory for {year}: {e}")
print("All downloads completed.")
def calculate_10yr_average_and_export():
rows, cols = 1200, 3600
total_cumulative = np.zeros((rows, cols), dtype=np.float32)
files = [f for f in os.listdir(BASE_DIR) if f.endswith(".dat.gz")]
for file in files:
# gsmap_gauge.YYYYMM.0.1d... から年と月を抽出
year_str = file.split(".")[1][0:4]
month_str = file.split(".")[1][4:6]
month_days = calendar.monthrange(int(year_str), int(month_str))[1]
filepath = os.path.join(BASE_DIR, file)
with gzip.open(filepath, 'rb') as f:
data = np.frombuffer(f.read(), dtype='<f4').reshape(2, rows, cols)[0]
data = np.where(data < 0, 0, data) # 欠損値処理
total_cumulative += data * 24 * month_days # mm/monthに換算して加算
# 10年分の総雨量を10で割り、年間平均降水量を算出
avg_annual = total_cumulative / (END_YEAR - START_YEAR + 1)
transform = from_origin(0.0, 60.0, 0.1, 0.1)
with rasterio.open(
OUTPUT_TIFF, 'w', driver='GTiff',
height=rows, width=cols, count=1, dtype=avg_annual.dtype,
crs='+proj=latlong +datum=WGS84', transform=transform, nodata=-9999
) as dst:
dst.write(avg_annual, 1)
print(f"Exported to {OUTPUT_TIFF}")
if __name__ == "__main__":
download_10years_data()
calculate_10yr_average_and_export()
Ａ part of OutputはBelowでございます。Very LongなOutputなので、AlmostOmitしています。
(Omit Above Here）
Downloading gsmap_gauge.202503.0.1d.monthly.v8.1000.0.dat.gz...
Downloading gsmap_gauge.202509.0.1d.monthly.v8.1000.0.dat.gz...
All downloads completed.
Exported to gsmap_10yr_avg_annual_precipitation.tif
これでPastTenYearsのRainFallのDataをGetしたうえで、10yearsのAnnualAvarageRainfallのGlobalなDistributionをProcess＆Culculatできました。
However、ThisCodeのままだと、DataがGoogleColabのTemporaryFolderにSaveされ、GoogleColabをCloseするためにMakeDeleteするので、TemporaryFolderにSaveされたFileをMyGoogleDriveのDirectedFolderにTransferしてDataを永続化しておきました。
それでは10YearsAverageAnnualRainfallをGlobalMapにDrawしてみます。CodeはNextです。
# ==============================================================================
# [Feature]
# Google Drive上のGeoTIFFファイル(10年平均の全球年間降水量)を読み込み、Matplotlibを用いて
# 地球全体の降水量分布をカラーマップとして図化します。欠損値を透過処理し、見やすい色域で描画します。
#
# [License & Disclaimer]
# Copyright (c) 2026 [Your Name/Organization] / Released under the MIT License.
# 【免責事項】本スクリプトによる可視化結果は衛星観測データに基づく推計値であり、実際の降水量や
# 水資源量を完全に保証するものではありません。本コードの利用により生じた直接的・間接的な損害について、
# 開発者は一切の責任を負いません。公開アプリ等へ組み込む際は十分な検証を行ってください。
# ==============================================================================
import rasterio
import numpy as np
import matplotlib.pyplot as plt
FILE_PATH = "/content/drive/MyDrive/Colab Notebooks/2026_JAXA_SateliteDataMeasere/2026_RainFallData/gsmap_10yr_avg_annual_precipitation.tif"
def plot_global_rainfall():
with rasterio.open(FILE_PATH) as src:
data = src.read(1)
# 欠損値(-9999や極端な負の値)をNaNに変換して描画対象から除外し、透過させる
data = np.where(data <= 0, np.nan, data)
# ラスターの座標系(Extent)を取得して軸に反映
extent = [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top]
plt.figure(figsize=(15, 7))
# 'Blues'カラーマップを使用し、上限(vmax)を3000mmに設定してコントラストを強調
im = plt.imshow(data, cmap='Blues', extent=extent, vmin=0, vmax=3000, origin='upper')
plt.colorbar(im, label='Average Annual Precipitation (mm/year)', fraction=0.02, pad=0.04)
plt.title("Global 10-Year Average Annual Precipitation (GSMaP)", fontsize=16)
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
# グラフの余白を自動調整して表示
plt.tight_layout()
plt.show()
if __name__ == "__main__":
plot_global_rainfall()
!pip install cartopy
# ==============================================================================
# [Feature]
# Google Drive上の全球年間降水量(GeoTIFF)を読み込み、経度を0〜360度から-180〜180度に
# シフトさせた上で、Cartopyを用いて白地図(海岸線・国境線)と重ね合わせて描画します。
#
# [License & Disclaimer]
# Copyright (c) 2026 [Your Name/Organization] / Released under the MIT License.
# 【免責事項】本スクリプトによる可視化結果は衛星観測データに基づく推計値であり、実際の降水量や
# 水資源量を完全に保証するものではありません。本コードの利用により生じた直接的・間接的な損害について、
# 開発者は一切の責任を負いません。公開アプリ等へ組み込む際は十分な検証を行ってください。
# ==============================================================================
import rasterio
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
FILE_PATH = "/content/drive/MyDrive/Colab Notebooks/2026_JAXA_SateliteDataMeasere/2026_RainFallData/gsmap_10yr_avg_annual_precipitation.tif"
def plot_global_rainfall_with_map():
with rasterio.open(FILE_PATH) as src:
data = src.read(1)
data = np.where(data <= 0, np.nan, data) # 欠損値を透過
# 経度(0〜360)を白地図(-180〜180)に合わせるため、3600ピクセルの半分(1800)を左右ローテート
data_shifted = np.roll(data, shift=1800, axis=1)
extent = [-180, 180, -60, 60]
# CartopyのPlateCarree(正距円筒図法)プロジェクションを指定
fig, ax = plt.subplots(figsize=(15, 7), subplot_kw={'projection': ccrs.PlateCarree()})
# 白地図（海岸線と国境線）を背景に追加
ax.add_feature(cfeature.COASTLINE, linewidth=0.8, edgecolor='black')
ax.add_feature(cfeature.BORDERS, linewidth=0.4, edgecolor='gray', linestyle=':')
# 降水量データを重ね合わせ
im = ax.imshow(data_shifted, cmap='Blues', extent=extent, vmin=0, vmax=3000,
origin='upper', transform=ccrs.PlateCarree())
plt.colorbar(im, ax=ax, label='Average Annual Precipitation (mm/year)', fraction=0.02, pad=0.04)
ax.set_title("Global 10-Year Average Annual Precipitation (GSMaP)", fontsize=16)
# 緯度・経度のグリッド線を追加
gl = ax.gridlines(draw_labels=True, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
gl.top_labels = False
gl.right_labels = False
plt.show()
if __name__ == "__main__":
plot_global_rainfall_with_map()
あと、VerticalAxisをLookすると、なぜかLatitudeがSouth、Northともに60DegreeでCutされています。このPointをSurveyすると、どうやらGSMaPのDataはArtificialSatelliteのSpec上North-SouthともにSurveyRangeがlessThan60Degreeまでのようです。ちょっと見にくいので、HigherLatitudeもMap上にDrawするように変更します。ただしNaturallyながら、HigherLatitudeのRainFallDataはNo-Existです。
SameDataをUseしてMatplotlibの別のHeatMapを使ってDrawしてみます。Colorによって見え方がかなり異なっているのがわかります。
- Blues（標準的な青色スケール）: 降水量といえば「水」。直感的に最も理解しやすく、ニュースや天気予報でもよく使われる王道の配色です。
- viridis（科学的標準スケール）: Matplotlibのデフォルトであり、視覚的な均一性が高く、色覚多様性（色盲など）にも配慮されたユニバーサルデザインの配色です。
- turbo（高コントラストスケール）: 虹色（jet）の欠点を改良した鮮やかな配色です。少雨地帯と多雨地帯のコントラストが極めて強く出るため、極端な気候の違いを強調したい場合に適しています。
MatPlotLibのOtherのHeatMapをTestするために、NextはAroundJapanのAnnualRainfallをDrawしてみます。TestはBelowのThreeKindsのHeatMapです。
- YlGnBu（イエロー・グリーン・ブルー）: 今回の土壌水分量でも使用した配色です。乾燥（黄）から湿潤（青）への変化が自然で、水資源の分析において最も違和感なく受け入れられやすい標準的なスケールです。
- magma（高コントラスト・ダーク系）: 黒〜紫〜オレンジ〜白へと変化します。背景が暗く沈むため、紀伊半島や四国南部など、降水量が極端に多い（オレンジ〜白になる）「局地的な豪雨エリア」をピンポイントで際立たせるのに非常に優れています。
- Spectral_r（赤〜黄〜青のダイバージェング）: 気候区分などでよく使われる配色です。少雨地帯が「赤〜オレンジ（暖色）」、多雨地帯が「青（寒色）」で表現されるため、直感的に「乾いた地域」と「濡れた地域」の境界を視覚化できます。
結論から申し上げますと、これは私たちの書いたコードの処理ミス（バグ）ではなく、JAXAの「GSMaP」データそのものに元から含まれている観測・補正アーティファクト（ノイズ）です。
衛星による降水観測データにこのような「四角い特異点」が現れるのには、主に以下の2つの科学的な理由があります。
地上観測データ（雨量計）の補正による歪み: GSMaPデータ（特に今回取得元のベースとなっているGSMaP_Gaugeというプロダクト）は、衛星の観測値を実際の地上の雨量計データとすり合わせて精度を上げています。しかし、中国の一部山岳地帯や北朝鮮などの「地上観測データが極端に少ない、あるいは国際的に共有されていない地域」では、1つの観測所のデータが数十km〜数百km四方の広範囲なグリッド（四角形）に強制的に適用されてしまい、隣のグリッドとの間に不自然な段差（ブロック状のノイズ）が生まれてしまうことがあります。
積雪・氷と降水の区別エラー: 衛星に搭載されているマイクロ波センサーは、地表に降り積もった雪や氷と、空から降っている最中の雨や雪を区別するのが非常に苦手です。冬場に雪に覆われる高緯度の山岳地帯などでは、地表の氷を「強い降水」として誤認してしまい、極端な異常値が出ることがあります。
逆に言えば、日本列島の降水分布が非常に滑らかで高精細に描画されているのは、日本全国に「アメダス」や気象レーダーといった世界トップクラスの緻密な地上観測網が張り巡らされており、衛星データの補正が完璧に機能している証拠でもあります。
I'm sure.
しかし、AroundJapanってもっとRainfallがあるのかと思っていましたが、Korean Peninsuraの北東あたりの方がRainがFallするということをあらためて知りました。
ちょっと、GSMaPから独自に計算した年間平均降雨量の精度を確認しておこうと思い、気象庁が発表している平年値と比べておくことにしました。GeminiにBelowCodeをWriteしてもらいました。
# ==============================================================================
# [Feature]
# 気象庁(JMA)のアメダス平年値と、JAXAのGSMaP衛星観測データ(10年平均)を比較し、
# 日本の主要5都市における年間降水量の精度検証(グラウンド・トゥルース)を行います。
# ==============================================================================
import rasterio
import numpy as np
import matplotlib.pyplot as plt
RAIN_FILE = "/content/drive/MyDrive/Colab Notebooks/2026_JAXA_SateliteDataMeasere/2026_RainFallData/gsmap_10yr_avg_annual_precipitation.tif"
# 気象庁の平年値(1991〜2020年)と各都市の代表緯度経度
city_data = {
"Sapporo": {"lat": 43.06, "lon": 141.33, "jma_rain": 1159.9},
"Tokyo": {"lat": 35.69, "lon": 139.75, "jma_rain": 1598.2},
"Niigata": {"lat": 37.90, "lon": 139.05, "jma_rain": 1806.9},
"Osaka": {"lat": 34.68, "lon": 135.52, "jma_rain": 1338.3},
"Fukuoka": {"lat": 33.58, "lon": 130.40, "jma_rain": 1612.3}
}
def plot_ground_truth_comparison():
cities = list(city_data.keys())
jma_values = [data["jma_rain"] for data in city_data.values()]
gsmap_values = []
# GSMaPデータから各都市のピクセル値を抽出
with rasterio.open(RAIN_FILE) as src:
coords = [(data["lon"], data["lat"]) for data in city_data.values()]
# src.sample はジェネレータを返すためループで取得
for val in src.sample(coords):
rain_val = val[0]
# 欠損値(-9999)の場合は0に丸める
gsmap_values.append(rain_val if rain_val != -9999 else 0.0)
# 棒グラフの描画設定
x = np.arange(len(cities))
width = 0.35
fig, ax = plt.subplots(figsize=(10, 6))
rects1 = ax.bar(x - width/2, jma_values, width, label='JMA (AMeDAS 30yr Normal)', color='#2ca02c')
rects2 = ax.bar(x + width/2, gsmap_values, width, label='JAXA (GSMaP 10yr Avg)', color='#1f77b4')
ax.set_ylabel('Annual Precipitation (mm/year)', fontsize=12)
ax.set_title('Ground Truth: Satellite Observation vs Ground Station', fontsize=14)
ax.set_xticks(x)
ax.set_xticklabels(cities, fontsize=12)
ax.legend(fontsize=11)
ax.grid(axis='y', linestyle='--', alpha=0.7)
# グラフ上に数値を表示
ax.bar_label(rects1, padding=3, fmt='%.0f')
ax.bar_label(rects2, padding=3, fmt='%.0f')
plt.tight_layout()
plt.show()
if __name__ == "__main__":
plot_ground_truth_comparison()
OutputはLikeBelowです。
気象庁の年間降水量平年値（1991〜2020年）は、私が事前に学習している気象庁公開の公式統計データ（「過去の気象データ検索」の平年値テーブル）から引用しました。
各管区気象台における公式な年間降水量の平年値は以下の通りです。
札幌: 1159.9 mm、東京: 1598.2 mm、新潟: 1806.9 mm、大阪: 1338.3 mm。福岡: 1612.3 mm
Excelでよく使うグラデーションがかかったBarGraphにして、少しDecorativeなGraphにTransformしてみます。MatPlotlibにはグラデーションの塗りつぶしがないそうで、Geminiがこんな工夫をしてくれました。
横方向にグラデーションをかけて円柱のような質感（ハイライトとシャドウ）を出し、平面のまま立体感を演出するアプローチですね！
Matplotlibの標準機能には「棒グラフのグラデーション塗り」がないため、1本の棒グラフを50分割の細いスライスに分け、端を暗く、中央やや左側を明るく（光が当たっているように）計算して描画するカスタム関数を実装しました。
# ==============================================================================
# [Feature]
# 気象庁(JMA)とJAXA(GSMaP)の降水量比較データを、横方向のグラデーションを
# 用いて円柱風の立体的な棒グラフとして描画します。
# ==============================================================================
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.colors as mcolors
def draw_cylinder_bars(ax, x_positions, heights, width, base_color, label):
"""棒を細いスライスに分割し、明るさを変えることで円柱風のグラデーションを描画します"""
n_slices = 50
slice_width = width / n_slices
rgb = np.array(mcolors.to_rgb(base_color))
# 凡例とデータラベル表示用に、透明(alpha=0)のベース棒グラフを先に描画
base_bars = ax.bar(x_positions, heights, width=width, color=base_color, label=label, alpha=0)
# グラデーションの描画
for x_pos, height in zip(x_positions, heights):
start_x = x_pos - width / 2
for i in range(n_slices):
# -1(左端) から 1(右端) までの相対位置
rel_pos = (i / (n_slices - 1)) * 2 - 1
# 左側(-0.2の位置)にハイライトが来るように明るさをカーブさせる
brightness = 1.0 - 0.4 * (rel_pos + 0.2)**2
color = np.clip(rgb * brightness, 0, 1)
ax.bar(start_x + i * slice_width, height, width=slice_width,
color=color, edgecolor='none', align='edge')
return base_bars
def plot_cylindrical_gradient_chart():
cities = ['Sapporo', 'Tokyo', 'Niigata', 'Osaka', 'Fukuoka']
jma_values = [1160, 1598, 1807, 1338, 1612]
gsmap_values = [1272, 1648, 1868, 1606, 2018]
x = np.arange(len(cities))
width = 0.35
fig, ax = plt.subplots(figsize=(10, 6))
# グラデーション円柱グラフの描画
bars1 = draw_cylinder_bars(ax, x - width/2, jma_values, width, '#2ca02c', 'JMA (AMeDAS 30yr Normal)')
bars2 = draw_cylinder_bars(ax, x + width/2, gsmap_values, width, '#1f77b4', 'JAXA (GSMaP 10yr Avg)')
# グラフの装飾
ax.set_ylabel('Annual Precipitation (mm/year)', fontsize=12)
ax.set_title('Ground Truth: Satellite Observation vs Ground Station (Cylindrical Style)', fontsize=14)
ax.set_xticks(x)
ax.set_xticklabels(cities, fontsize=12)
# 凡例とグリッドの設定
legend = ax.legend(fontsize=11)
for handle in legend.legend_handles:
handle.set_alpha(1.0) # 凡例の色を不透明に戻す
ax.grid(axis='y', linestyle='--', alpha=0.7)
# 数値ラベルの追加
ax.bar_label(bars1, padding=3, fmt='%.0f')
ax.bar_label(bars2, padding=3, fmt='%.0f')
plt.tight_layout()
plt.show()
if __name__ == "__main__":
plot_cylindrical_gradient_chart()
Nextに、DataについてCosiderationでございます。
全体としてJAXAのDataがJMAのDataよりも大きいですね。さらに大阪、福岡では300～400㎜の誤差があるのが、とても気になります。技術的、あるいは論理的な理由はありますか？
衛星観測（面）と地上観測（点）で生じるこの誤差は、リモートセンシング分野において非常にクラシックかつ重要な議論の的です。全体的なJAXAデータの過大評価、および特に大阪や福岡で300〜400mmの乖離が出ている現象には、明確な技術的・気象学的な理由が4つ存在します。
観測期間のズレ（気候変動の影響）: 今回比較した気象庁のアメダスデータは「過去30年間（1991〜2020年）」の平年値ですが、GSMaPデータは「直近10年間（2016〜2025年）」の平均です。近年激増している線状降水帯やゲリラ豪雨、台風の大型化により、直近10年の方が純粋に降水量が多くなっているという気象トレンドが直接反映されています。
「点」と「面」の空間解像度の違い: アメダスは直径わずか20cmの雨量計の「点」の観測ですが、GSMaPは約10km四方の「面」の平均値です。例えば夏の局地的なゲリラ豪雨がアメダスの観測所を数キロ逸れて降った場合、アメダスは「降水ゼロ」を記録しますが、10kmメッシュの衛星はその雨雲をしっかり拾ってカウントします。
沿岸部と複雑な地形によるマイクロ波の誤認: 大阪（大阪湾と六甲・生駒山系）や福岡（玄界灘と背振山地）は、海と山が極めて近接した地形です。土壌水分マップの時にも直面しましたが、海陸が混在する沿岸ピクセルや山岳の急斜面では、マイクロ波の反射特性が複雑になり、アルゴリズムが氷の粒などを「より強い雨」として過大評価（オーバーエスティメート）しやすい弱点を持っています。
強風による雨量計の「アンダーキャッチ（捕捉率低下）」: 台風や冬型の強い季節風が吹く時、地上の雨量計は雨粒や雪が風で横に流されてしまい、実際に降った量よりも10〜20%ほど少なく計測してしまう物理的な欠点があります。宇宙から雨雲全体を捉える衛星にはこの風による漏れがありません。
総じて、衛星データは地形の影響で「過大評価」しやすく、地上データは強風などの影響で「過小評価」しやすいという相反する特性を持っています。これらの特性を理解した上で両者をブレンドしていくのが、現代の気象データ解析の基本スタンスです。
End
コメント
コメントを投稿