#! /usr/bin/env python3
#
#
# Purpose:
# Test manipulation of .grb2 files from NOAA WAVEWATCH III hindcast repo
#
# Notes:
# 1. pygrib is not compatible with Windows.
# 2. GRIB = General Regularly distributed Information in Binary form
# GRIB is a binary format of gridded data.
# 3. WAVEWATCH III: https://polar.ncep.noaa.gov/waves/hindcasts/nopp-phase2.php
# 3. pygrib docs: https://jswhit.github.io/pygrib/api.html#example-usage
#
# Author(s):
# Caleb Grant
#
# Revision History
# Date Reason
# ----------------------------
# 2021-09-23 Created. CG.
# ===========================================GRIB
Manipulating gridded datasets
# Standard libraries
import datetime
import os
import numpy as np
# Third party packages
import pygrib
# Make output of plotting commands display inline with notebook
%matplotlib inline--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) Cell In[2], line 6 3 import datetime 5 # Third party packages ----> 6 import pygrib 7 import numpy as np 9 # Make output of plotting commands display inline with notebook ModuleNotFoundError: No module named 'pygrib'
Test scenario
data_dir = os.path.join("data", "200901")
# U-component of input wind, V-component of input wind (m/s).
# Vector described by these components shows direction in which whind blows
wind_grb_file = "multi_reanal.ecg_10m.wind.200901.grb2"
wind_grb_obj = pygrib.open(os.path.join(data_dir, wind_grb_file))
# Significant height of combined wind waves and swell (m)
hs_grb_file = "multi_reanal.ecg_10m.hs.200901.grb2"
hs_grb_obj = pygrib.open(os.path.join(data_dir, hs_grb_file))
# Primary wave mean period (s)
tp_grb_file = "multi_reanal.ecg_10m.tp.200901.grb2"
tp_grb_obj = pygrib.open(os.path.join(data_dir, tp_grb_file))
# Primary wave direction (degrees true, i.e. 0 deg = coming from North, 90 deg = coming from East)
dp_grb_file = "multi_reanal.ecg_10m.dp.200901.grb2"
dp_grb_obj = pygrib.open(os.path.join(data_dir, dp_grb_file))Object methods
print(dir(dp_grb_obj))['__call__', '__class__', '__delattr__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__iter__', '__le__', '__lt__', '__ne__', '__new__', '__next__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '_advance', 'close', 'closed', 'has_multi_field_msgs', 'message', 'messagenumber', 'messages', 'name', 'read', 'readline', 'rewind', 'seek', 'select', 'tell']GRIB Object keys
print(sorted(dp_grb_obj[1].keys()))['GRIBEditionNumber', 'NV', 'Ni', 'Nj', 'PLPresent', 'PVPresent', 'alternativeRowScanning', 'analDate', 'angleDivisor', 'angleMultiplier', 'angleSubdivisions', 'average', 'backgroundProcess', 'basicAngleOfTheInitialProductionDomain', 'binaryScaleFactor', 'bitMapIndicator', 'bitmapPresent', 'bitsPerValue', 'bottomLevel', 'centre', 'centreDescription', 'cfName', 'cfNameECMF', 'cfVarName', 'cfVarNameECMF', 'changeDecimalPrecision', 'codedValues', 'dataDate', 'dataRepresentationTemplateNumber', 'dataTime', 'day', 'decimalPrecision', 'decimalScaleFactor', 'deleteCalendarId', 'deleteLocalDefinition', 'deletePV', 'discipline', 'distinctLatitudes', 'distinctLongitudes', 'editionNumber', 'endStep', 'forecastTime', 'g2grid', 'genVertHeightCoords', 'generatingProcessIdentifier', 'getNumberOfValues', 'globalDomain', 'grib2LocalSectionPresent', 'grib2divider', 'gridDefinitionDescription', 'gridDefinitionTemplateNumber', 'gridDescriptionSectionPresent', 'gridType', 'hour', 'hoursAfterDataCutoff', 'iDirectionIncrement', 'iDirectionIncrementGiven', 'iDirectionIncrementInDegrees', 'iScansNegatively', 'iScansPositively', 'identifier', 'ieeeFloats', 'ifsParam', 'ijDirectionIncrementGiven', 'indicatorOfUnitOfTimeRange', 'interpretationOfNumberOfPoints', 'isConstant', 'isHindcast', 'is_aerosol', 'is_aerosol_optical', 'is_chemical', 'is_chemical_distfn', 'is_efas', 'is_uerra', 'jDirectionIncrement', 'jDirectionIncrementGiven', 'jDirectionIncrementInDegrees', 'jPointsAreConsecutive', 'jScansPositively', 'julianDay', 'kurtosis', 'latLonValues', 'latitudeOfFirstGridPoint', 'latitudeOfFirstGridPointInDegrees', 'latitudeOfLastGridPoint', 'latitudeOfLastGridPointInDegrees', 'latitudes', 'lengthOfHeaders', 'level', 'localTablesVersion', 'longitudeOfFirstGridPoint', 'longitudeOfFirstGridPointInDegrees', 'longitudeOfLastGridPoint', 'longitudeOfLastGridPointInDegrees', 'longitudes', 'mAngleMultiplier', 'mBasicAngle', 'masterDir', 'maximum', 'md5Headers', 'md5Section1', 'md5Section3', 'md5Section4', 'md5Section5', 'md5Section6', 'md5Section7', 'minimum', 'minute', 'minutesAfterDataCutoff', 'missingValue', 'modelName', 'month', 'name', 'nameECMF', 'nameOfFirstFixedSurface', 'nameOfSecondFixedSurface', 'neitherPresent', 'numberOfDataPoints', 'numberOfMissing', 'numberOfOctetsForNumberOfPoints', 'numberOfSection', 'numberOfSection', 'numberOfSection', 'numberOfSection', 'numberOfSection', 'numberOfSection', 'numberOfValues', 'offsetValuesBy', 'optimizeScaleFactor', 'packingType', 'paramId', 'paramIdECMF', 'parameterCategory', 'parameterName', 'parameterNumber', 'parameterUnits', 'pressureUnits', 'productDefinitionTemplateNumber', 'productType', 'productionStatusOfProcessedData', 'radius', 'referenceValue', 'referenceValueError', 'resolutionAndComponentFlags', 'resolutionAndComponentFlags1', 'resolutionAndComponentFlags2', 'resolutionAndComponentFlags6', 'resolutionAndComponentFlags7', 'resolutionAndComponentFlags8', 'scaleFactorOfEarthMajorAxis', 'scaleFactorOfEarthMinorAxis', 'scaleFactorOfFirstFixedSurface', 'scaleFactorOfRadiusOfSphericalEarth', 'scaleFactorOfSecondFixedSurface', 'scaleValuesBy', 'scaledValueOfEarthMajorAxis', 'scaledValueOfEarthMinorAxis', 'scaledValueOfFirstFixedSurface', 'scaledValueOfRadiusOfSphericalEarth', 'scaledValueOfSecondFixedSurface', 'scanningMode', 'scanningMode5', 'scanningMode6', 'scanningMode7', 'scanningMode8', 'second', 'section0Length', 'section1Length', 'section3Length', 'section4Length', 'section5Length', 'section6Length', 'section7Length', 'section8Length', 'sectionNumber', 'sectionNumber', 'sectionNumber', 'sectionNumber', 'sectionNumber', 'sectionNumber', 'sectionNumber', 'selectStepTemplateInstant', 'selectStepTemplateInterval', 'setBitsPerValue', 'setCalendarId', 'shapeOfTheEarth', 'shortName', 'shortNameECMF', 'significanceOfReferenceTime', 'skewness', 'sourceOfGridDefinition', 'standardDeviation', 'startStep', 'stepRange', 'stepType', 'stepTypeInternal', 'stepUnits', 'subCentre', 'subdivisionsOfBasicAngle', 'tablesVersion', 'tablesVersionLatest', 'targetCompressionRatio', 'tempPressureUnits', 'topLevel', 'totalLength', 'typeOfCompressionUsed', 'typeOfFirstFixedSurface', 'typeOfGeneratingProcess', 'typeOfLevel', 'typeOfOriginalFieldValues', 'typeOfProcessedData', 'typeOfSecondFixedSurface', 'units', 'unitsECMF', 'unitsOfFirstFixedSurface', 'unitsOfSecondFixedSurface', 'uvRelativeToGrid', 'validDate', 'validityDate', 'validityTime', 'values', 'year']Test value arrays
values = np.array([grb.values for grb in hs_grb_obj])
print(values.shape, values.min(), values.max())(249, 331, 301) 0.0 9999.0print(dir(values))['T', '__abs__', '__add__', '__and__', '__array__', '__array_finalize__', '__array_function__', '__array_interface__', '__array_prepare__', '__array_priority__', '__array_struct__', '__array_ufunc__', '__array_wrap__', '__bool__', '__class__', '__complex__', '__contains__', '__copy__', '__deepcopy__', '__delattr__', '__delitem__', '__dir__', '__divmod__', '__doc__', '__eq__', '__float__', '__floordiv__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__iadd__', '__iand__', '__ifloordiv__', '__ilshift__', '__imatmul__', '__imod__', '__imul__', '__index__', '__init__', '__init_subclass__', '__int__', '__invert__', '__ior__', '__ipow__', '__irshift__', '__isub__', '__iter__', '__itruediv__', '__ixor__', '__le__', '__len__', '__lshift__', '__lt__', '__matmul__', '__mod__', '__mul__', '__ne__', '__neg__', '__new__', '__or__', '__pos__', '__pow__', '__radd__', '__rand__', '__rdivmod__', '__reduce__', '__reduce_ex__', '__repr__', '__rfloordiv__', '__rlshift__', '__rmatmul__', '__rmod__', '__rmul__', '__ror__', '__rpow__', '__rrshift__', '__rshift__', '__rsub__', '__rtruediv__', '__rxor__', '__setattr__', '__setitem__', '__setstate__', '__sizeof__', '__str__', '__sub__', '__subclasshook__', '__truediv__', '__xor__', 'all', 'any', 'argmax', 'argmin', 'argpartition', 'argsort', 'astype', 'base', 'byteswap', 'choose', 'clip', 'compress', 'conj', 'conjugate', 'copy', 'ctypes', 'cumprod', 'cumsum', 'data', 'diagonal', 'dot', 'dtype', 'dump', 'dumps', 'fill', 'flags', 'flat', 'flatten', 'getfield', 'imag', 'item', 'itemset', 'itemsize', 'max', 'mean', 'min', 'nbytes', 'ndim', 'newbyteorder', 'nonzero', 'partition', 'prod', 'ptp', 'put', 'ravel', 'real', 'repeat', 'reshape', 'resize', 'round', 'searchsorted', 'setfield', 'setflags', 'shape', 'size', 'sort', 'squeeze', 'std', 'strides', 'sum', 'swapaxes', 'take', 'tobytes', 'tofile', 'tolist', 'tostring', 'trace', 'transpose', 'var', 'view']# sets, rows per set, columns per row
print(values.shape)
# would equate to
sets = len(values)
rows = len(values[0])
cols = len(values[0][0])
print(sets, rows, cols)(249, 331, 301)
249 331 301Testing lat/longs
for grb in dp_grb_obj[1:6]:
lats, lons = grb.latlons()
print("min/max lat ang lon: ", lats.min(), lats.max(), lons.min(), lons.max())min/max lat ang lon: -0.00010999999961602835 55.0 260.0 310.0000000000057
min/max lat ang lon: -0.00010999999961602835 55.0 260.0 310.0000000000057
min/max lat ang lon: -0.00010999999961602835 55.0 260.0 310.0000000000057
min/max lat ang lon: -0.00010999999961602835 55.0 260.0 310.0000000000057
min/max lat ang lon: -0.00010999999961602835 55.0 260.0 310.0000000000057wind_grb_file = "multi_reanal.ecg_10m.wind.200901.grb2"
wind_grb_obj = pygrib.open(os.path.join(data_dir, wind_grb_file))
# Reqind the iterator
wind_grb_obj.rewind()
# values = np.array([grb.values for grb in wind_grb_obj])
# np.savetxt("test_wind.csv", values, delimiter=",")
data_dict = {}
# Each iteration captures point in time across entire grid
for grb in wind_grb_obj:
data_dict.update({grb.validDate: {grb.parameterName: np.array(grb.values)}})print(data_dict[datetime.datetime(2009, 1, 1, 0, 0)]){'v-component of wind': array([[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
...,
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.]])}hs_grb_obj.rewind()
print(hs_grb_obj[1].parameterName)
date_selected = datetime.datetime(2009, 1, 1, 0, 0)
values = []
for grb in hs_grb_obj:
if grb.validDate == date_selected and grb.parameterName == "Significant height of combined wind waves and swell":
values.append(grb.values)
values = np.array(values)
print(values.shape, values.min(), values.max())
lats, longs = grb.latlons()
print(lats.min(), lats.max(), longs.min(), longs.max())Significant height of combined wind waves and swell
(1, 331, 301) 0.0 9999.0
-0.00010999999961602835 55.0 260.0 310.0000000000057"""
fig = plt.figure(figsize=(16,35))
# https://matplotlib.org/basemap/users/ortho.html
m = Basemap(projection='lcc', lon_0=-80, lat_0=30, resolution='l', width=5.e6, height=5.e6)
x, y = m(longs, lats)
# for vals in range(1, 2):
# print(values[vals])
ax = plt.plot(values[0])
m.drawcoastlines()
cs = m.contourf(x, y, values[0], np.linspace(230, 300, 41), cmap=plt.cm.jet, extend='both')
t = plt.title('test')
""""\nfig = plt.figure(figsize=(16,35))\n# https://matplotlib.org/basemap/users/ortho.html\nm = Basemap(projection='lcc', lon_0=-80, lat_0=30, resolution='l', width=5.e6, height=5.e6)\nx, y = m(longs, lats)\n\n# for vals in range(1, 2):\n# print(values[vals])\n\nax = plt.plot(values[0])\nm.drawcoastlines()\ncs = m.contourf(x, y, values[0], np.linspace(230, 300, 41), cmap=plt.cm.jet, extend='both')\nt = plt.title('test')\n"hs_grb_obj.rewind()
for grb in hs_grb_obj[1:2]:
print(grb.validDate, grb.parameterName, grb.latlons())2009-01-01 00:00:00 Significant height of combined wind waves and swell (array([[ 5.5000000e+01, 5.5000000e+01, 5.5000000e+01, ...,
5.5000000e+01, 5.5000000e+01, 5.5000000e+01],
[ 5.4833333e+01, 5.4833333e+01, 5.4833333e+01, ...,
5.4833333e+01, 5.4833333e+01, 5.4833333e+01],
[ 5.4666666e+01, 5.4666666e+01, 5.4666666e+01, ...,
5.4666666e+01, 5.4666666e+01, 5.4666666e+01],
...,
[ 3.3322400e-01, 3.3322400e-01, 3.3322400e-01, ...,
3.3322400e-01, 3.3322400e-01, 3.3322400e-01],
[ 1.6655700e-01, 1.6655700e-01, 1.6655700e-01, ...,
1.6655700e-01, 1.6655700e-01, 1.6655700e-01],
[-1.1000000e-04, -1.1000000e-04, -1.1000000e-04, ...,
-1.1000000e-04, -1.1000000e-04, -1.1000000e-04]]), array([[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
...,
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ]]))print((datetime.datetime(2009, 1, 1, 0, 0) + datetime.timedelta(hours=744)).strftime("%Y-%m-%d"))2009-02-01Better iteration method:
# hs_grb_obj.rewind()
hs_grb_obj.seek(0)
print(hs_grb_obj.tell())
print(hs_grb_obj.read(1))0
[1:Significant height of combined wind waves and swell:m (instant):regular_ll:surface:level 1:fcst time 0 hrs:from 200901010000]hs_grb_obj.seek(0)
print(hs_grb_obj.tell())
for grb in hs_grb_obj:
grb
print(hs_grb_obj.tell())0
249hs_grb_obj.seek(0)
grb = hs_grb_obj.read(1)[0]
print("1")
print(grb.values)
print("2")
grb = hs_grb_obj.select(name="Significant height of combined wind waves and swell")[0]
print(grb.values)1
[[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
...
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]]
2
[[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
...
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]]hs_grb_obj.select(name="Significant height of combined wind waves and swell")[1:6][2:Significant height of combined wind waves and swell:m (instant):regular_ll:surface:level 1:fcst time 3 hrs:from 200901010000,
3:Significant height of combined wind waves and swell:m (instant):regular_ll:surface:level 1:fcst time 6 hrs:from 200901010000,
4:Significant height of combined wind waves and swell:m (instant):regular_ll:surface:level 1:fcst time 9 hrs:from 200901010000,
5:Significant height of combined wind waves and swell:m (instant):regular_ll:surface:level 1:fcst time 12 hrs:from 200901010000,
6:Significant height of combined wind waves and swell:m (instant):regular_ll:surface:level 1:fcst time 15 hrs:from 200901010000]Find the first grib message with a matching name:
grb = hs_grb_obj.select(name="Significant height of combined wind waves and swell")[0]Extract the data values using the values key (grb.keys() will return a list of the available keys):
# The data is returned as a numpy array, or if missing values or a bitmap
# are present, a numpy masked array. Reduced lat/lon or gaussian grid
# data is automatically expanded to a regular grid. Details of the internal
# representation of the grib data (such as the scanning mode) are handled
# automatically.
grb.valuesmasked_array(
data=[[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
...,
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --]],
mask=[[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
...,
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True]],
fill_value=9999.0)grb.values.shape, grb.values.min(), grb.values.max()((331, 301), 0.0, 6.5600000000000005)Get the latitudes and longitudes of the grid:
lats, lons = grb.latlons()
lats.shape, lats.min(), lats.max(), lons.shape, lons.min(), lons.max()((331, 301),
-0.00010999999961602835,
55.0,
(331, 301),
260.0,
310.0000000000057)Extract data and get lat/lon values for a subset over North America:
data, lats, lons = grb.data(lat1=20, lat2=70, lon1=220, lon2=320)
data.shape, lats.min(), lats.max(), lons.min(), lons.max()((210, 301), 20.166597000000415, 55.0, 260.0, 310.0000000000057)Close the grib file
hs_grb_obj.close()Testing docstrings
grbs = wind_grb_obj
grbs.seek(0)
u_grb = grbs.select(name="U component of wind")
v_grb = grbs.select(name="V component of wind")
print(len(u_grb), len(v_grb))249 249print(u_grb[0].values.shape, v_grb[0].values.shape)(331, 301) (331, 301)grbs = wind_grb_obj
grbs.seek(0)
grb = grbs.select(name="U component of wind")[0]grb.has_key("gaulats")Falsepygrib.gribmessage object
print("grb.messagenumber: ", grb.messagenumber)
print("grb.fcstimeunits: ", grb.fcstimeunits)
print("grb.analDate: ", grb.analDate)
print("grb.validDate: ", grb.validDate)grb.messagenumber: 1
grb.fcstimeunits: hrs
grb.analDate: 2009-01-01 00:00:00
grb.validDate: 2009-01-01 00:00:00# If the default values of lat1,lat2,lon1,lon2 are None, which means the entire grid is returned
grb.data(lat1=None, lat2=None, lon1=None, lon2=None)(masked_array(
data=[[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
...,
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --]],
mask=[[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
...,
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True]],
fill_value=9999.0),
array([[ 5.5000000e+01, 5.5000000e+01, 5.5000000e+01, ...,
5.5000000e+01, 5.5000000e+01, 5.5000000e+01],
[ 5.4833333e+01, 5.4833333e+01, 5.4833333e+01, ...,
5.4833333e+01, 5.4833333e+01, 5.4833333e+01],
[ 5.4666666e+01, 5.4666666e+01, 5.4666666e+01, ...,
5.4666666e+01, 5.4666666e+01, 5.4666666e+01],
...,
[ 3.3322400e-01, 3.3322400e-01, 3.3322400e-01, ...,
3.3322400e-01, 3.3322400e-01, 3.3322400e-01],
[ 1.6655700e-01, 1.6655700e-01, 1.6655700e-01, ...,
1.6655700e-01, 1.6655700e-01, 1.6655700e-01],
[-1.1000000e-04, -1.1000000e-04, -1.1000000e-04, ...,
-1.1000000e-04, -1.1000000e-04, -1.1000000e-04]]),
array([[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
...,
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ],
[260. , 260.16666667, 260.33333333, ..., 309.66666667,
309.83333333, 310. ]]))grbs = wind_grb_obj
grbs.seek(0)
u_grb = grbs.select(name="U component of wind")
v_grb = grbs.select(name="V component of wind")
grbs = u_grbfor grb in grbs:
print("grb.messagenumber: ", grb.messagenumber)
print("grb.fcstimeunits: ", grb.fcstimeunits)
print("grb.validDate: ", grb.validDate)
print("type(grb.values): ", type(grb.values))
print(grb.values)
breakgrb.messagenumber: 1
grb.fcstimeunits: hrs
grb.validDate: 2009-01-01 00:00:00
type(grb.values): <class 'numpy.ma.core.MaskedArray'>
[[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
...
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]
[-- -- -- ... -- -- --]]Masked arrays
# Masked arrays are arrays that may have missing or invalid entries
# https://numpy.org/doc/stable/reference/maskedarray.generic.html
# When an element of the mask is False, the corresponding element of the associated array is valid and is said to be unmasked.
# When an element of the mask is True, the corresponding element of the associated array is said to be masked (invalid).
# The package ensures that masked entries are not used in computations.
import numpy.ma as maMasked array example
x = np.array([1, 2, 3, -1, 5])
print(x.mean())
# Mark the fourth element as invalid
mx = ma.masked_array(x, mask=[0, 0, 0, 1, 0])
print(mx.mean())2.0
2.75Usage
grbs = wind_grb_obj
grbs.seek(0)
u_grb = grbs.select(name="U component of wind")
v_grb = grbs.select(name="V component of wind")
grb = u_grb[0]# Test if input is instance of masked array
ma.isMaskedArray(grb.values)True# Check if all elements evaluate to True
# If false, invalid values exist
ma.all(grb.values)Falseprint(ma.count(grb.values)) # non-maked elements
print(ma.count_masked(grb.values)) # masked elements30096
69535ma.getdata(grb.values)array([[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
...,
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.]])ma.ravel(grb.values) # Returns a 1D version of self, as a view.masked_array(data=[--, --, --, ..., --, --, --],
mask=[ True, True, True, ..., True, True, True],
fill_value=9999.0)ma.ravel(grb.values).fill_value9999.0ma.filled(grb.values)array([[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
...,
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.],
[9999., 9999., 9999., ..., 9999., 9999., 9999.]])ma.MaskedArray.torecords(grb.values)array([[(9999., True), (9999., True), (9999., True), ...,
(9999., True), (9999., True), (9999., True)],
[(9999., True), (9999., True), (9999., True), ...,
(9999., True), (9999., True), (9999., True)],
[(9999., True), (9999., True), (9999., True), ...,
(9999., True), (9999., True), (9999., True)],
...,
[(9999., True), (9999., True), (9999., True), ...,
(9999., True), (9999., True), (9999., True)],
[(9999., True), (9999., True), (9999., True), ...,
(9999., True), (9999., True), (9999., True)],
[(9999., True), (9999., True), (9999., True), ...,
(9999., True), (9999., True), (9999., True)]],
dtype=[('_data', '<f8'), ('_mask', '?')])wind_grb_obj.name # name of object
wind_grb_obj.messages # count of messages
wind_grb_obj.messagenumber # Current message index0