Retrieve and process monthly rainfall maps retrieved from GSMaP¶

Notebook prepared by Ben Maathuis. ITC-University of Twente, Enschede. The Netherlands

An important source of rainfall (accumulated for different aggregation periods - Algorithm Version 6) as well as climatology and standardized precipitation indexdata, available from April 2000 can be retrieved from the JAXA Global Rainfall Watch portal, available at: https://sharaku.eorc.jaxa.jp/GSMaP/. In this notebook use we are going to download monthly precipitation and Standarized Precipitation Index data for a given timestamp, so the methodology can be used to create a longer time series.

For additional background information see: "Global Satellite Mapping of Precipitation (GSMaP) Products in the GPM Era" (https://link.springer.com/chapter/10.1007/978-3-030-24568-9_20). Also available at the GSMap SFTP site (see folder \docs)

The link to the monthly averaged rain rate data: /climate/gnrt6/monthly/YYYY/gsmap_gnrt6.YYYYMM.0.1d.monthly.dat.gz

The link to the monthly averaged rain rate data: /climate/gnrt6/monthly/YYYY/gsmap_gnrt6.YYYYMM.0.25d.monthly.spi01.dat.gz where:

  • YYYY: 4-digit year;
  • MM: 2-digit month;

If not done already, complete the (free) registration to get access to the data archive, see: https://sharaku.eorc.jaxa.jp/GSMaP/registration.html, your Username and Password are required in the code field below

In [1]:
#paramiko is a Python package for working with SSH, SFTP, and SCP-like file transfers. 
#It's commonly used to automate connections to Linux servers, HPC clusters, and remote file systems.

#!pip install paramiko
In [2]:
import paramiko
import numpy as np
#import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import ilwis
import os
import struct
import gzip
import shutil
import warnings
warnings.filterwarnings("ignore")
In [3]:
#data folder
work_dir = os.getcwd()+r'/GSMap'
work_dir = os.path.normpath(work_dir)

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir) 

if os.path.isdir(work_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(work_dir)

current dir is: d:\jupyter\notebook_scripts\Special\GSMaP_rainfall
current working directory is: d:\jupyter\notebook_scripts\Special\GSMaP_rainfall\GSMap
Folder exists
In [4]:
#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)

d:\jupyter\notebook_scripts\Special\GSMaP_rainfall\GSMap
In [5]:
#provide date for monthly PCP layer - filename convention = gsmap_gnrt6.YYYYMM.0.1d.monthly.dat.gz
pyear = 2010
pmonth = 12
In [6]:
#prepare string for data retrieval given the date provided
yy = str(pyear)
mm = str(pmonth).zfill(2) #month always with 2 digits 


filename = (
    "/climate/gnrt6/monthly/"
    + str(pyear)
    + "/gsmap_gnrt6."
    + yy
    + mm
    + ".0.1d.monthly.dat.gz"
)

print(filename)

/climate/gnrt6/monthly/2010/gsmap_gnrt6.201012.0.1d.monthly.dat.gz
In [7]:
#Provide your GSMap registration credentials
USERNAME = "your username"
PASSWORD = "your password"
In [8]:
#download the data - here monthly PCP

HOST = "hokusai.eorc.jaxa.jp"
PORT = 22


REMOTE_FILE = filename

LOCAL_FILE = filename.split("/")[-1] #only filename is used

ssh = paramiko.SSHClient()
ssh.set_missing_host_key_policy(paramiko.AutoAddPolicy())

try:
    ssh.connect(
        hostname=HOST,
        port=PORT,
        username=USERNAME,
        password=PASSWORD
    )

    sftp = ssh.open_sftp()

    sftp.get(REMOTE_FILE, work_dir+'/'+LOCAL_FILE)

    print(f"Downloaded: {LOCAL_FILE}")

finally:
    sftp.close()
    ssh.close()

Downloaded: gsmap_gnrt6.201012.0.1d.monthly.dat.gz
In [9]:
#unzip the downloaded file
gz_file = work_dir+'/'+LOCAL_FILE
out_file = work_dir+'/'+"gsmap_"+yy+mm+"_pcp.dat"

with gzip.open(gz_file, "rb") as f_in:
    with open(out_file, "wb") as f_out:
        shutil.copyfileobj(f_in, f_out)

print("Decompressed to:", out_file)

Decompressed to: d:\jupyter\notebook_scripts\Special\GSMaP_rainfall\GSMap/gsmap_201012_pcp.dat

GSMap monthly rainfall: All binary files are produced in little-endian byte order platform, and archived with compressed using “gzip”. In each monthly file, there are two global fields: monthly averaged rain rate; and numbers of valid pixel (≧ 0 mm) per month. The former unit is [mm/hr] and the missing value is -999.9. Multiplying of both layers gives the monthly total precipitation [mm/month].

In [10]:
# read GSMap monthly file content into memory

cols = 3600
rows = 1200
bands = 2 
bytesperpixel = 4 # float

f = open(out_file, 'rb')
buffer = f.read(rows * cols * bands * bytesperpixel)
it = struct.iter_unpack("<f", buffer)
data_pcp = np.array([el[0] for el in it])
In [11]:
data_pcp.shape
Out[11]:
(8640000,)
In [12]:
data_3d = data_pcp.reshape(2, 1200, 3600)
data_3d.shape
Out[12]:
(2, 1200, 3600)
In [13]:
#display the numpy array created in the previous steps using matplotlib 
fig = plt.figure(figsize =(6, 6))
plt.imshow(data_3d[0], interpolation='none', vmin=0, vmax=1.5, cmap='jet')
Out[13]:
<matplotlib.image.AxesImage at 0x27f673dea50>
No description has been provided for this image
In [14]:
# Make a float copy so NaN can be stored
layer1 = data_3d[0].astype(float)
layer2 = data_3d[1].astype(float)

# Replace the missing value
layer1[layer1 == -999] = np.nan

# Multiply
result = layer1 * layer2
In [15]:
#Split into western and eastern halves and concatonate
western = result[:, :1800]
eastern = result[:, 1800:]

data_shifted = np.concatenate((eastern, western), axis=1)
In [16]:
lon = np.arange(-180, 180, 0.1)
lat = np.arange(-60, 60, 0.1)

print(len(lon))  # 3600
print(len(lat))  # 1200

3600
1200
In [17]:
#display the Monthly PCP map

fig = plt.figure(figsize=(12, 4))
ax = plt.axes(projection=ccrs.PlateCarree())

im = ax.imshow(
    data_shifted,
    extent=[lon.min(), lon.max(), lat.min(), lat.max()],
    transform=ccrs.PlateCarree(),
    #aspect="auto",
    vmin= 0,
    vmax= 500,
    cmap="jet"
)

# Add map features
ax.coastlines(resolution="110m",  linewidth=0.3, color="black")
ax.add_feature(cfeature.BORDERS,  linewidth=0.3)

# Borders in black
ax.add_feature(
    cfeature.BORDERS,
    edgecolor="black",
    linewidth=0.3
)


# Gridlines
gl = ax.gridlines(
    draw_labels=True,
    linewidth=0.5,
    alpha=0.5
)
gl.top_labels = False
gl.right_labels = False

cbar = plt.colorbar(im, ax=ax, pad=0.02, shrink=0.65)
cbar.set_label("Monthly PCP (mm)")

ax.set_title("GSMaP Monthly PCP for "+str(pyear)+str(mm))

plt.show()
No description has been provided for this image

Export as ILWIS map¶

In [18]:
# envelope =  min X, max Y, max X, min Y,  e.g. -180 60 180 -60 (units - degree)
# gridsize =  no columns, no rows = 3600 / 1200
grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=-180 60 180 -60 , gridsize=3600 1200, cornerofcorners=yes')
In [19]:
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 10000.0, 0))
pcpNew = ilwis.RasterCoverage()
pcpNew.setSize(ilwis.Size(3600, 1200,1))
pcpNew.setGeoReference(grf_LL)
pcpNew.setDataDef(dfNum)
In [20]:
# 1 D array
p1 = np.array([data_shifted]).flatten()
print(p1.shape)
print(p1)

(4320000,)
[-283971.60693359 -284971.50695801       0.         ... -435956.41064453
 -437956.21069336 -432956.71057129]
In [21]:
pcpNew.array2raster(p1)
print(pcpNew.size())

Size(3600, 1200, 1)
In [22]:
# Store monthly PCP map
pcpNew1 = ilwis.do('mapcalc', 'iff(@1<=0, ?, @1)',pcpNew) 
pcpNew1.store('GSMaP_'+str(pyear)+str(mm)+'_pcp.mpr')

Retrieve and process monthly SPI3 maps retrieved from GSMaP¶

The 3-month Standardized Precipitation Index (SPI-3) is a meteorological drought indicator that measures precipitation anomalies over a 90-day period. It is primarily used to monitor short-to-medium-term moisture conditions, acting as an early proxy for agricultural impacts, reduced soil moisture, and changes in smaller creek flows. The index calculates how observed precipitation deviates from a long-term historical average for the same specific 3-month window. Because it is normalized, SPI values represent standard deviations from the mean:

  • Positive values (1.0 or higher) indicate wetter-than-normal conditions.
  • Near normal conditions fall between ‒0.99 and 0.99.
  • Negative values (‒1.0 or lower) indicate a precipitation deficit. For instance, a value from ‒1.0 to ‒1.49 signifies "moderately dry," while values below ‒2.0 indicate "extremely dry" or severe drought.

SPI Archived data: /climate/gnrt6/SPI/YYYY/gsmap_gnrt6.YYYYMM.0.25d.monthly.spi0Z.dat.gz

where,

  • YYYY: 4-digit year;
  • MM: 2-digit month;
  • Z: 1-digit number (1-3). Note here the 3 month SPI is used

We are going to retrieve the SPI data for the same timestep as the precipitation data processed before.

In [23]:
#prepare string for data retrieval given the date provided
yy = str(pyear)
mm = str(pmonth).zfill(2) #month always with 2 digits 


filename = (
    "/climate/gnrt6/SPI/"
    + str(pyear)
    + "/gsmap_gnrt6."
    + yy
    + mm
    + ".0.25d.monthly.spi03.dat.gz"
)

print(filename)

/climate/gnrt6/SPI/2010/gsmap_gnrt6.201012.0.25d.monthly.spi03.dat.gz
In [24]:
#download the data - here 3-month SPI

HOST = "hokusai.eorc.jaxa.jp"
PORT = 22

REMOTE_FILE = filename

LOCAL_FILE = filename.split("/")[-1] #only filename is used

ssh = paramiko.SSHClient()
ssh.set_missing_host_key_policy(paramiko.AutoAddPolicy())

try:
    ssh.connect(
        hostname=HOST,
        port=PORT,
        username=USERNAME,
        password=PASSWORD
    )

    sftp = ssh.open_sftp()

    sftp.get(REMOTE_FILE, work_dir+'/'+LOCAL_FILE)

    print(f"Downloaded: {LOCAL_FILE}")

finally:
    sftp.close()
    ssh.close()

Downloaded: gsmap_gnrt6.201012.0.25d.monthly.spi03.dat.gz
In [25]:
gz_file = work_dir+'/'+LOCAL_FILE
out_file = work_dir+'/'+"gsmap_"+yy+mm+"_spi03.dat"

with gzip.open(gz_file, "rb") as f_in:
    with open(out_file, "wb") as f_out:
        shutil.copyfileobj(f_in, f_out)

print("Decompressed to:", out_file)

Decompressed to: d:\jupyter\notebook_scripts\Special\GSMaP_rainfall\GSMap/gsmap_201012_spi03.dat
In [26]:
# read SPI3 file  content into memory

cols = 1440
rows = 480
bands = 1 
bytesperpixel = 4 # float

#file = work_dir +'/PRCP_CU_GAUGE_V1.0GLB_0.50deg.lnx.20110103.RT'
f = open(out_file, 'rb')
buffer = f.read(rows * cols * bands * bytesperpixel)
it = struct.iter_unpack("<f", buffer)
data_spi3 = np.array([el[0] for el in it])
In [27]:
data_spi3.shape
Out[27]:
(691200,)
In [28]:
data_3d = data_spi3.reshape(480, 1440)
data_3d.shape
Out[28]:
(480, 1440)
In [29]:
#quick display the numpy array created in the previous steps using matplotlib 
fig = plt.figure(figsize =(6, 6))
plt.imshow(data_3d, interpolation='none')
Out[29]:
<matplotlib.image.AxesImage at 0x27f0008c410>
No description has been provided for this image
In [30]:
#Split into western and eastern halves and concatonate
western = data_3d[:, :720]
eastern = data_3d[:, 720:]

data_shifted = np.concatenate((eastern, western), axis=1)
In [31]:
lon = np.arange(-180, 180, 0.25)
lat = np.arange(-60, 60, 0.25)

print(len(lon))  # 1440
print(len(lat))  # 480

1440
480
In [32]:
#display the SPI3 map
data_shifted[data_shifted == -999] = np.nan
fig = plt.figure(figsize=(12, 4))
ax = plt.axes(projection=ccrs.PlateCarree())

im = ax.imshow(
    data_shifted,
    extent=[lon.min(), lon.max(), lat.min(), lat.max()],
    transform=ccrs.PlateCarree(),
    aspect="auto",
    vmin= -4,
    vmax= 4,
    cmap="jet_r"
)

# Add map features
ax.coastlines(resolution="110m",  linewidth=0.3, color="black")
ax.add_feature(cfeature.BORDERS,  linewidth=0.3)

# Borders in black
ax.add_feature(
    cfeature.BORDERS,
    edgecolor="black",
    linewidth=0.3
)


# Gridlines
gl = ax.gridlines(
    draw_labels=True,
    linewidth=0.5,
    alpha=0.5
)
gl.top_labels = False
gl.right_labels = False

cbar = plt.colorbar(im, ax=ax, pad=0.02)
cbar.set_label("SPI3")

ax.set_title("GSMaP SPI3 for "+str(pyear)+str(mm))

plt.show()
No description has been provided for this image

SPI Index values:

  • Extremely Wet: 2.0+
  • Very Wet: 1.5 to 1.99
  • Moderately Wet: 1.0 to 1.49
  • Near Normal: -0.99 to 0.99
  • Moderately Dry: -1.0 to -1.49
  • Severely Dry: -1.5 to -1.99
  • Extremely Dry: -2.0 or less
In [33]:
# envelope =  min X, max Y, max X, min Y,  e.g. -180 60 180 -60 (units - degree)
# gridsize =  no columns, no rows = 1440 / 480
grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=-180 60 180 -60 , gridsize=1440 480, cornerofcorners=yes')
In [34]:
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-10.0, 10.0, 0))
spiNew = ilwis.RasterCoverage()
spiNew.setSize(ilwis.Size(1440, 480,1))
spiNew.setGeoReference(grf_LL)
spiNew.setDataDef(dfNum)
In [35]:
# 1 D array
p1 = np.array([data_shifted]).flatten()
print(p1.shape)
print(p1)

(691200,)
[nan nan nan ... nan nan nan]
In [36]:
spiNew.array2raster(p1)
print(spiNew.size())

Size(1440, 480, 1)
In [37]:
# Store spi3 map
spiNew1 = ilwis.do('mapcalc', 'iff(@1==-999, ?, @1)',spiNew) 
spiNew1.store('GSMaP_'+str(pyear)+str(mm)+'_spi3.mpr')
In [ ]: