Analyse de l'agence météorologique Les précipitations (format GRIB2) sont gérées en Python sous la forme d'un tableau Numpy

Je veux lire l'analyse des précipitations en Python

Le fichier bin du maillage de 1 km pluie analysée (2006-) fourni par l'Agence météorologique / Centre de soutien météorologique aux entreprises est au format GRIB2. Cependant, à partir d'avril 2020, il ne peut plus être géré par GDAL ou pygrib car il contient ses propres extensions. Il existe également une méthode de vidage binaire utilisant wgrib2 qui prend en charge l'analyse des précipitations ou la convertit en NetCDF et la lit, mais si vous souhaitez extraire uniquement la valeur numérique, il est plus pratique de la lire directement.

Analyse de la structure du fichier GRIB2 des précipitations

Les données se composent de ** 8 ** (sections 0 à 8, dont la section 2 est omise), dont la valeur de niveau de l'entier représentant les précipitations dans la section 7 est la longueur de la course comprimée à l'ouest. 2560 éléments sont stockés du nord à l'est [^ 1]. Les précipitations sont converties en convertissant cette valeur de niveau en fonction d'une table de recherche avec la ** valeur représentative des données ** (10 fois les précipitations représentant la grille de 1 km) définie dans la section 5, et en la multipliant par 1/10. Tu peux l'avoir.

Chaque section de données

Les spécifications des données binaires sont les Informations techniques sur les matériels de distribution (édition météorologique) n ° 238 et n ° 193 édition révisée Il est décrit dans /238.pdf). Les modèles utilisés dans chaque section sont spécifiés dans International Weather Report Ceremony / Separate Volume.

Chaque section a une section de longueur fixe et une section de longueur variable, et la longueur de section (nombre d'octets) et le numéro de section sont stockés au début de la section. Pour lire une section particulière, commencez la lecture en ajoutant les longueurs de la section précédente. Dans le format de données d'analyse des précipitations, les sections 0, 1, 3, 4 et 6 sont définies comme des longueurs fixes de 16, 21, 72, 82 et 6 octets dans l'ordre, et les sections 5 et 7 sont le nombre de niveaux et de précipitations. Il est variable selon la distribution. Cependant, dans la section 5, l'analyse actuelle des précipitations n'a pas changé à un maximum de 98 niveaux, elle est donc en fait de 213 octets (= 17 + 2 * 98).

Valeur représentative des données

Le tableau ci-dessous montre la correspondance entre les valeurs de niveau et les valeurs représentatives des données définies dans la section 5.

En plus de cela, la section 7 inclut le niveau 0 comme «valeur manquante». Lors du décodage en tant que valeur entière, il est pratique d'attribuer une valeur appropriée, donc dans le code ci-dessous, elle est définie sur "-10".

Codage de la longueur de tirage

La méthode de codage de longueur d'exécution utilisée dans la section 7 exprime le nombre de répétitions par la somme des nombres, et [Explication de la méthode de codage de longueur d'exécution](http://www.jmbsc.or.jp/jp/jp Développez en fonction de la formule (/online/joho-sample/Run-Length_Encoding.pdf).

Référence: Toyoda, E., 2012: modèles GRIB2 pour les données de précipitations analysées par radar / pluviomètre JMA: PDT 4.50008 et DRT 5.200

Avant l'expansion, il s'agit d'une chaîne d'octets d'entiers non signés (0 à 255). La valeur de niveau de 0 à la valeur la plus élevée (<= 98) ou moins qui apparaît est la valeur de niveau telle quelle et la valeur de 0 à 255 correspondant au nombre de répétitions de la valeur de niveau précédent. Si ce dernier est continu, la somme des nombres ajoutés comme nombre de chiffres élevés dans l'ordre est le nombre de répétitions.

Décoder

Voici la quantité de précipitations contenue dans l'échantillon du fichier d'analyse des précipitations téléchargé à partir du Meteorological Business Support Center. Code à lire comme ndarray de Numpy (Nord-Sud 3360 x Est-Ouest 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()

Dessiner avec Matplotlib et Basemap

Pour le tableau bidimensionnel «pluie» des précipitations obtenues, générez la latitude et la longitude et dessinez-les avec Matplotlib. Ici, considérant que les données sont une valeur représentative dans la grille, les coordonnées sont divisées en 80 parties égales de la grille de maillage cubique (80 parties égales de 1 degré dans le sens de la longitude et 40 minutes (= 1 / 1,5 degré) dans le sens de la latitude). ) Est défini pour être le centre (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()

Si vous développez un fichier, ce sera 65,6 Mo dans np.float64. Si vous renvoyez la valeur telle qu'elle est 10 fois, ce sera environ 16,4 Mo de 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]: Dans le 72e octet de la section 4, 0x00 est spécifié comme mode de balayage, qui est` les premières données à l'extrémité nord-ouest de la plage, empruntant la représentation de format.txt attachée à l'analyse des précipitations. Oui, si vous vous alignez d'ouest en est à la même latitude et que vous avancez du nombre de points de grille le long de la ligne de latitude, vous vous déplacerez vers les données à l'extrémité ouest d'une grille au sud, et les données à la même latitude seront à nouveau alignées par le nombre de points de grille le long de la ligne de latitude.