Analyse der Meteorologischen Agentur Niederschlag (GRIB2-Format) wird in Python als Numpy-Array behandelt

Ich möchte den Analyse-Niederschlag in Python lesen

Die vom Meteorological Agency / Meteorological Business Support Center bereitgestellte Bin-Datei des 1 km-Netzes analysierter Niederschlag (2006-) liegt im GRIB2-Format vor. Ab April 2020 kann es jedoch nicht von GDAL oder pygrib verarbeitet werden, da es eigene Erweiterungen enthält. Es gibt auch eine Methode zum binären Dumping mit wgrib2, die den Analyse-Niederschlag unterstützt oder ihn in NetCDF konvertiert und liest. Wenn Sie jedoch nur den numerischen Wert extrahieren möchten, ist es bequemer, ihn direkt zu lesen.

Analyse Niederschlag GRIB2-Dateistruktur

Die Daten bestehen aus ** 8 ** (Abschnitte 0 bis 8, von denen Abschnitt 2 weggelassen ist), von denen der Pegelwert der Ganzzahl, die den Niederschlag in Abschnitt 7 darstellt, nach Westen komprimiert ist. 2560 Elemente werden von Nord nach Ost gespeichert [^ 1]. Der Niederschlag wird konvertiert, indem dieser Pegelwert basierend auf einer Nachschlagetabelle mit dem in Abschnitt 5 definierten ** Datenrepräsentativwert ** (10-facher Niederschlag, der das 1 km-Gitter darstellt) konvertiert und mit 1/10 multipliziert wird. Sie können es durch bekommen

Jeder Datenabschnitt

Die Spezifikationen der Binärdaten sind die Technischen Informationen zu Verteilungsmaterialien (Meteorologische Ausgabe) Nr. 238 und Nr. 193 überarbeitete Ausgabe. Es ist in /238.pdf) beschrieben. Die in den einzelnen Abschnitten verwendeten Vorlagen sind in Internationale Wetterberichtszeremonie / Separater Band angegeben.

Jeder Abschnitt hat einen Abschnitt fester Länge und einen Abschnitt variabler Länge, und die Abschnittslänge (Anzahl der Bytes) und die Abschnittsnummer werden am Anfang des Abschnitts gespeichert. Um einen bestimmten Abschnitt zu lesen, beginnen Sie mit dem Lesen aus der Summe der Längen des vorherigen Abschnitts. Im Analyse-Niederschlagsdatenformat sind die Abschnitte 0, 1, 3, 4 und 6 als feste Längen von 16, 21, 72, 82 und 6 Bytes in der Reihenfolge definiert, und die Abschnitte 5 und 7 sind die Anzahl der Ebenen und der Niederschlag. Sie ist je nach Verteilung variabel. Abschnitt 5 umfasst jedoch tatsächlich 213 Byte (= 17 + 2 * 98), da sich der aktuelle Analyseniederschlag auf dem Höchststand von 98 nicht geändert hat.

Datenrepräsentativer Wert

Die folgende Tabelle zeigt die Entsprechung zwischen den in Abschnitt 5 definierten Pegelwerten und datenrepräsentativen Werten.

Darüber hinaus enthält Abschnitt 7 Stufe 0 als "fehlenden Wert". Bei der Dekodierung als ganzzahliger Wert ist es zweckmäßig, einen geeigneten Wert zuzuweisen, sodass dieser im folgenden Code auf "-10" gesetzt wird.

Lauflängencodierung

Die in Abschnitt 7 verwendete Lauflängencodierungsmethode drückt die Anzahl der Wiederholungen durch die Summe der Zahlen und Erläuterung der Lauflängencodierungsmethode aus Erweitern Sie basierend auf der Formel (/online/joho-sample/Run-Length_Encoding.pdf).

Referenz: Toyoda, E., 2012: GRIB2-Vorlagen für JMA-Radar- / Regenmesser-analysierte Niederschlagsdaten: PDT 4.50008 und DRT 5.200

Vor der Erweiterung handelt es sich um eine Bytefolge von Ganzzahlen ohne Vorzeichen (0 bis 255). Der angezeigte Pegelwert von 0 bis zum höchsten Wert (<= 98) oder weniger ist der aktuelle Pegelwert und der Wert von 0 bis 255 entspricht der Anzahl der Wiederholungen des vorherigen Pegelwerts. Wenn letzteres stetig ist, ist die Summe der Zahlen, die als Anzahl der in der Reihenfolge angehobenen Ziffern addiert werden, die Anzahl der Wiederholungen.

Dekodieren

Das Folgende ist der Niederschlag, der in der Stichprobe der Analyse-Niederschlagsdatei enthalten ist, die vom Meteorological Business Support Center heruntergeladen wurde. Code als Numpys Ndarray (Nord-Süd 3360 x Ost-West 2560).

from itertools import repeat
import struct
import numpy as np

def set_table(section5):
    max_level = struct.unpack_from('>H', section5, 15)[0]
    table = (
        -10, # define representative of level 0　(Missing Value)
        *struct.unpack_from('>'+str(max_level)+'H', section5, 18)
    )
    return np.array(table, dtype=np.int16)

def decode_runlength(code, hi_level):
    for raw in code:
        if raw <= hi_level:
            level = raw
            pwr = 0
            yield level
        else:
            length = (0xFF - hi_level)**pwr * (raw - (hi_level + 1))
            pwr += 1
            yield from repeat(level, length)

def load_jmara_grib2(file):
    with open(file, 'rb') as f:
        binary = f.read()

    len_ = {'sec0':16, 'sec1':21, 'sec3':72, 'sec4':82, 'sec6':6}
    
    end4 = len_['sec0'] + len_['sec1'] + len_['sec3'] + len_['sec4'] - 1
    len_['sec5'] = struct.unpack_from('>I', binary, end4+1)[0]
    section5 = binary[end4:(end4+len_['sec5']+1)]

    end6 = end4 + len_['sec5'] + len_['sec6']
    len_['sec7'] = struct.unpack_from('>I', binary, end6+1)[0]
    section7 = binary[end6:(end6+len_['sec7']+1)]
    
    highest_level = struct.unpack_from('>H', section5, 13)[0]
    level_table = set_table(section5)
    decoded = np.fromiter(
        decode_runlength(section7[6:], highest_level), dtype=np.int16
    ).reshape((3360, 2560))

    # convert level to representative
    return level_table[decoded]

file = rb'./SRF_GPV_Ggis1km_Prr60lv_ANAL_20180707/Z__C_RJTD_20180707000000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'
rain = load_jmara_grib2(file) / 10

subroutine load_jmara_grib2(filepath, decoded_2d)
    implicit none
    character(*), intent(in) :: filepath
    integer, parameter :: len_sec0 = 16, len_sec1 = 21, len_sec2 = 0
    integer, parameter :: len_sec3 = 72, len_sec4 = 82, len_sec6 = 6
    integer :: len_sec5, len_sec7, pos, filesize

    integer, parameter :: x_n = 2560, y_n = 3360
    integer(1), allocatable :: binary(:)
    integer(2), allocatable :: table(:)
    integer(2) :: highest_level, decoded(x_n*y_n)
    integer(2), intent(out) :: decoded_2d(x_n, y_n)

    inquire(file=filepath, size=filesize)
    allocate(binary(filesize))
    open(100, file=trim(filepath), status='old', access='direct', form='unformatted', recl=filesize)
    read(100, rec=1) binary
    close(100)

    pos = len_sec0 + len_sec1 + len_sec2 + len_sec3 + len_sec4
    len_sec5 = pack_int32(binary(pos+1), binary(pos+2), binary(pos+3), binary(pos+4))
    highest_level = pack_int16(binary(pos+13), binary(pos+14))
    call set_table(binary(pos+1 : pos+len_sec5), table)

    pos = pos + len_sec5 + len_sec6
    len_sec7 = pack_int32(binary(pos+1), binary(pos+2), binary(pos+3), binary(pos+4))
    call decode_runlength(binary(pos+6 : pos+len_sec7), highest_level, decoded)

    where (decoded == 0)
        decoded = -10   ! define representative of level 0 (Missing Value)
    elsewhere
        decoded = table(decoded)
    endwhere

    decoded_2d = reshape(decoded, (/ x_n, y_n /))

contains

    subroutine set_table(section5, table)
        implicit none
        integer(1), intent(in) :: section5(:)
        integer(2), allocatable, intent(out) :: table(:)
        integer(2) :: max_level, level

        max_level = pack_int16(section5(15), section5(16))
        allocate(table(max_level))
        do level = 1, max_level
            table(level) = pack_int16(section5(16 + 2*level), section5(17 + 2*level))
        end do
    end subroutine

    subroutine decode_runlength(code, hi_level, decoded)
        implicit none
        integer(1), intent(in) :: code(:)
        integer(2), intent(in) :: hi_level
        integer(2), intent(out) :: decoded(0:)
        integer :: i, pwr, offset, length
        integer(2) :: raw, level
        
        offset = 0
        do i = 1, size(code)
            raw = get_int16(code(i))
            if(raw <= hi_level) then
                level = raw
                pwr = 0
                decoded(offset) = level
                offset = offset + 1
            else
                length = (255 - hi_level)**pwr * (raw - (hi_level + 1))
                pwr = pwr + 1
                decoded(offset : offset+length-1) = level
                offset = offset + length
            end if
        end do
    end subroutine

    integer(2) function pack_int16(arg1, arg2)
        implicit none
        integer(1) :: arg1, arg2

        pack_int16 = ishft( iand(int(B'11111111',2), int(arg1, 2)), 8 )     &
                   +        iand(int(B'11111111',2), int(arg2, 2))
    end function

    integer(4) function pack_int32(arg1, arg2, arg3, arg4)
        implicit none
        integer(1) :: arg1, arg2, arg3, arg4

        pack_int32 = ishft( iand(int(B'11111111',4), int(arg1, 4)), 8*3)    &
                   + ishft( iand(int(B'11111111',4), int(arg2, 4)), 8*2)    &
                   + ishft( iand(int(B'11111111',4), int(arg3, 4)), 8  )    &
                   +        iand(int(B'11111111',4), int(arg4, 4))
    end function

    integer(2) function get_int16(arg1)
        implicit none
        integer(1) :: arg1

        get_int16 = iand(int(B'11111111',2), int(arg1, 2))
    end function
end subroutine

gfortran -shared jmara.f90 -Ofast -march=native -fopt-info -o ./jmara.dll

import ctypes
import numpy as np

def load_jmara_grib2(file):
    lib = np.ctypeslib.load_library('jmara.dll', '.')

    lib.load_jmara_grib2_.argtypes = [
        ctypes.c_char_p,
        np.ctypeslib.ndpointer(dtype=np.int16),
        ctypes.c_size_t     # hidden argument for character length
        ]

    lib.load_jmara_grib2_.restype = ctypes.c_void_p
    representative = np.zeros((3360, 2560), dtype=np.int16, order='C')

    lib.load_jmara_grib2_(
        ctypes.c_char_p(file),
        representative,
        ctypes.c_size_t(len(file))
        )
    
    return representative

file = rb'./SRF_GPV_Ggis1km_Prr60lv_ANAL_20180707/Z__C_RJTD_20180707000000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'
rain = load_jmara_grib2(file) / 10

import matplotlib.pyplot as plt

plt.figure(figsize=(12, 8))
cmap=plt.cm.jet
cmap.set_under('white')
plt.imshow(rain, aspect=8/12, vmin=0.4, cmap=cmap)
plt.colorbar()
plt.show()

Zeichnen Sie mit Matplotlib und Grundkarte

Generieren Sie für das zweidimensionale Array "Regen" des erhaltenen Niederschlags den Breiten- und Längengrad und zeichnen Sie ihn mit Matplotlib. Wenn man bedenkt, dass die Daten ein repräsentativer Wert im Gitter sind, werden die Koordinaten in 80 gleiche Teile des kubischen Maschengitters unterteilt (80 gleiche Teile von 1 Grad in Längengradrichtung und 40 Minuten (= 1 / 1,5 Grad) in Breitengradrichtung). ) Ist auf die Mitte eingestellt (1/2).

conda install -c conda-forge basemap basemap-data-hires

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

lon = np.linspace(118, 150, 2560, endpoint=False) + 1/80 / 2
lat = np.linspace(48, 20, 3360, endpoint=False) - 1/80/1.5 / 2
X, Y = np.meshgrid(lon, lat)

fig, ax = plt.subplots(figsize=(10,8))
ax.set_title('JMA Radar/Raingauge-Analyzed Precipitation\n2018-07-07T00:00:00+00:00')
m = Basemap(llcrnrlat=20, urcrnrlat=48, 
            llcrnrlon=118, urcrnrlon=150,
            resolution='l', ax=ax)
m.drawmeridians(np.arange(0, 360, 5), labels=[1, 0, 0, 1], linewidth=0)
m.drawparallels(np.arange(-90, 90, 5), labels=[1, 0, 1, 0], linewidth=0)
m.drawcoastlines(color='black')

levels = [0, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80]
colors = ['whitesmoke', 'gainsboro', 'darkgray', 'cyan', 'deepskyblue', 
          'blue', 'lime', 'yellow', 'orange', 'red']

im = m.contourf(X, Y, rain, levels, colors=colors, extend='both')
im.cmap.set_under('white')
im.cmap.set_over('magenta')
cb = m.colorbar(im, 'right', ticks=levels, size='2.5%')
cb.set_label('mm/hr')

plt.show()

Wenn Sie eine Datei erweitern, beträgt diese in np.float64 65,6 MB. Wenn Sie den Wert 10 Mal zurückgeben, beträgt er ungefähr 16,4 MB np.int16.

import datetime

def half_hour_range(start, end):
    for n in range(int((end - start).total_seconds() / 3600)):
        yield start + datetime.timedelta(hours=n)
        yield start + datetime.timedelta(hours=n, minutes=30)

def jmara_path(time, parent_dir='.'):
    path = parent_dir.rstrip('/').rstrip('\\') + (
        time.strftime('/%Y/DATA/%Y/%m/%d/')
        + 'Z__C_RJTD_'
        + time.strftime('%Y%m%d%H%M%S')
        + '_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'
    )
    return path

start = datetime.datetime(2018, 7, 1, 0, 0, 0)
end = datetime.datetime(2018, 9, 1, 0, 0, 0)
parent_dir = 'G:'

filepaths = [jmara_path(t, parent_dir) for t in half_hour_range(start, end)]

['G:/2018/DATA/2018/07/01/Z__C_RJTD_20180701000000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/07/01/Z__C_RJTD_20180701003000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/07/01/Z__C_RJTD_20180701010000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/07/01/Z__C_RJTD_20180701013000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/07/01/Z__C_RJTD_20180701020000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 (...)
 'G:/2018/DATA/2018/08/31/Z__C_RJTD_20180831213000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/08/31/Z__C_RJTD_20180831220000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/08/31/Z__C_RJTD_20180831223000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/08/31/Z__C_RJTD_20180831230000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin',
 'G:/2018/DATA/2018/08/31/Z__C_RJTD_20180831233000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin']

[^ 1]: Im 72. Oktett von Abschnitt 4 wird "0x00" als Scanmodus angegeben. Dies sind die ersten Daten am nordwestlichen Ende des Bereichs, wobei die dem Formatniederschlag beigefügte Darstellung von format.txt übernommen wird. Ja, wenn Sie von West nach Ost auf demselben Breitengrad ausgerichtet sind und um die Anzahl der Gitterpunkte entlang der Breitengradlinie vorrücken, werden Sie zu den Daten am westlichen Ende eines Gitters nach Süden verschoben, und die Daten auf demselben Breitengrad werden erneut durch die Anzahl der Gitterpunkte entlang der Breitengradlinie ausgerichtet. "