search

UMR-CNRM / EPyGrAM

Enhanced Python for Graphics and Analysis of Meteorological fields
Other
14 stars 15 forks source link

Problem when reading and rewriting a grib file #19

Open vincentchabot opened 10 months ago

vincentchabot commented 10 months ago

Hello, I have a problem when opening a grib file, reading a field and rewriting it. It seems that when writing the file it has some encoding problem. 


import os 
import shutil 
import epygram 
epygram.init_env()

print(epygram.__version__)
os.system("printenv")

fname =  "/cnrm/proc/chabotv/PanguWeather/NO_SAVE/input_data/grib/ARPEGE_2023100612_MER_GLOB025.grib"
ftmp = "test_input.grib"
# Copying the input file 
shutil.copy(fname,ftmp)
resAearo = epygram.formats.resource(filename=ftmp,openmode='a')
MSL = resAearo.readfield({
  'shortName': 'prmsl',
})
msl = MSL.getdata()
# Checking that data are correct
print(msl)
# Rewriting the file 
resAearo.writefield(MSL)
resAearo.close()
print("Writing has been done. Trying to read.")
resAearo = epygram.formats.resource(filename=ftmp,openmode='r')
MSL = resAearo.readfield({
  'shortName': 'prmsl',
})
msl = MSL.getdata()
print(msl)

Version of epygram 1.4.17

Environment variables MAGPLUS_DEBUG=off LESSOPEN=| /usr/bin/lesspipe %s FER_WEB_BROWSER=firefox CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt INPUTRC=/etc/inputrc MAIL=/var/spool/mail/chabotv FER_DATA=. /usr/share/ferret-vis/data /usr/share/ferret-vis/jnls/go /usr/share/ferret-vis/jnls/examples /usr/share/ferret-vis/jnls/contrib USER=chabotv NCARG_ROOT=/usr SSH_CLIENT=137.129.183.45 50448 22 LOCAL_DEFINITION_TEMPLATES=/usr/share/emos/gribtemplates/ HOSTNAME=sxgmap3 MARS_LSM_PATH=/usr/share/emos/interpol/ XDG_SESSION_TYPE=tty EPYGRAM_DIR=/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM GIT_SSL_CAINFO=/usr/local/etc/proxy1.pem SHLVL=2 LD_LIBRARY_PATH=/opt/compat-ubuntu/lib BROWSER=firefox MOTD_SHOWN=pam HOME=/home/chabotv OLDPWD=/home/chabotv XCDP_DEFAULT_USER=vincent SSH_TTY=/dev/pts/15 MV_EMOS_TABLE2_DIR=/usr/share/emos/gribtables/ FER_DESCR=. /usr/share/ferret-vis/descr DBUS_SESSION_BUS_ADDRESS=unix:path=/run/user/19171/bus FER_PALETTE=. /usr/share/ferret-vis/ppl/palettes ECMWF_LOCAL_TABLE_PATH=/usr/share/emos/gribtables/ LOGNAME=chabotv MAGPLUS_HOME=/usr/lib/python3/dist-packages/ecmwflibs _=/usr/bin/python3.10 MTOOLDIR=/cnrm/proc/chabotv/NO_SAVE XDG_SESSION_CLASS=user ECCODES_SAMPLES_PATH=/lib/share/eccodes/samples TERM=xterm-256color XDG_SESSION_ID=12701 FER_GRIDS=. /usr/share/ferret-vis/grids MAGPLUS_INFO=on HISTCONTROL=ignoreboth PATH=/home/common/sync/vortex/vortex-olive/bin:/home/common/sync/vortex/vortex-olive/site/arpifs_listings/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/opt/puppetlabs/bin:/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/apptools:/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/apptools FER_DIR=/usr/share/ferret-vis REQUESTS_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt METVIEW_EXTRA_GRIB_DEFINITION_PATH=/usr/share/eccodes/definitions LM_LICENSE_FILE=7496@flexa.cnrm.meteo.fr:7496@flexg.cnrm.meteo.fr:7496@flex.cnrm.meteo.fr XDG_RUNTIME_DIR=/run/user/19171 DISPLAY=sxgmap3:14.0 HISTSIZE=1000 LANG=C.UTF-8 LS_COLORS=rs=0:di=01;34:ln=01;36:mh=00:pi=40;33:so=01;35:do=01;35:bd=40;33;01:cd=40;33;01:or=40;31;01:mi=00:su=37;41:sg=30;43:ca=30;41:tw=30;42:ow=34;42:st=37;44:ex=01;32:*.tar=01;31:*.tgz=01;31:*.arc=01;31:*.arj=01;31:*.taz=01;31:*.lha=01;31:*.lz4=01;31:*.lzh=01;31:*.lzma=01;31:*.tlz=01;31:*.txz=01;31:*.tzo=01;31:*.t7z=01;31:*.zip=01;31:*.z=01;31:*.dz=01;31:*.gz=01;31:*.lrz=01;31:*.lz=01;31:*.lzo=01;31:*.xz=01;31:*.zst=01;31:*.tzst=01;31:*.bz2=01;31:*.bz=01;31:*.tbz=01;31:*.tbz2=01;31:*.tz=01;31:*.deb=01;31:*.rpm=01;31:*.jar=01;31:*.war=01;31:*.ear=01;31:*.sar=01;31:*.rar=01;31:*.alz=01;31:*.ace=01;31:*.zoo=01;31:*.cpio=01;31:*.7z=01;31:*.rz=01;31:*.cab=01;31:*.wim=01;31:*.swm=01;31:*.dwm=01;31:*.esd=01;31:*.jpg=01;35:*.jpeg=01;35:*.mjpg=01;35:*.mjpeg=01;35:*.gif=01;35:*.bmp=01;35:*.pbm=01;35:*.pgm=01;35:*.ppm=01;35:*.tga=01;35:*.xbm=01;35:*.xpm=01;35:*.tif=01;35:*.tiff=01;35:*.png=01;35:*.svg=01;35:*.svgz=01;35:*.mng=01;35:*.pcx=01;35:*.mov=01;35:*.mpg=01;35:*.mpeg=01;35:*.m2v=01;35:*.mkv=01;35:*.webm=01;35:*.webp=01;35:*.ogm=01;35:*.mp4=01;35:*.m4v=01;35:*.mp4v=01;35:*.vob=01;35:*.qt=01;35:*.nuv=01;35:*.wmv=01;35:*.asf=01;35:*.rm=01;35:*.rmvb=01;35:*.flc=01;35:*.avi=01;35:*.fli=01;35:*.flv=01;35:*.gl=01;35:*.dl=01;35:*.xcf=01;35:*.xwd=01;35:*.yuv=01;35:*.cgm=01;35:*.emf=01;35:*.ogv=01;35:*.ogx=01;35:*.aac=00;36:*.au=00;36:*.flac=00;36:*.m4a=00;36:*.mid=00;36:*.midi=00;36:*.mka=00;36:*.mp3=00;36:*.mpc=00;36:*.ogg=00;36:*.ra=00;36:*.wav=00;36:*.oga=00;36:*.opus=00;36:*.spx=00;36:*.xspf=00;36: BUFR_TABLES=/usr/share/emos/bufrtables/ FER_GO=. /usr/share/ferret-vis/jnls/go /usr/share/ferret-vis/jnls/examples /usr/share/ferret-vis/jnls/contrib SHELL=/bin/bash PROFILE_LOADED=1 LESSCLOSE=/usr/bin/lesspipe %s %s FER_DSETS=/usr/share/ferret-vis FONTCONFIG_FILE=/tmp/tmpw84ulpjpecmwflibs.xml PWD=/home/chabotv/tmp PROJ_LIB=/usr/lib/python3/dist-packages/ecmwflibs/share/proj SSH_CONNECTION=137.129.183.45 50448 137.129.166.118 22 XDG_DATA_DIRS=/usr/share/gnome:/usr/local/share:/usr/share:/var/lib/snapd/desktop LFI_HNDL_SPEC=:1 PYTHONPATH=/home/common/sync/vortex/vortex-olive:/home/common/sync/vortex/vortex-olive/src:/home/common/sync/vortex/vortex-olive/site::/home/chabotv/monorepo4ai/libs:/home/chabotv/monorepo4ai/projects/pangu_weather:/home/common/sync/vortex/vortex-olive:/home/common/sync/vortex/vortex-olive/src:/home/common/sync/vortex/vortex-olive/site::/home/chabotv/envar OMP_NUM_THREADS=1 FER_FONTS=/usr/share/ferret-vis/ppl/fonts MANPATH=:/opt/puppetlabs/puppet/share/man DR_HOOK_NOT_MPI=1 XCDP_DEFAULT_PASS=chabot 
[[ 99754.5234375  99754.5234375  99754.5234375 ...  99754.5234375
   99754.5234375  99754.5234375]
 [ 99950.5234375  99950.5234375  99950.5234375 ...  99950.5234375
   99950.5234375  99950.5234375]
 [100150.5234375 100150.5234375 100150.5234375 ... 100150.5234375
  100150.5234375 100150.5234375]
 ...
 [102594.5234375 102594.5234375 102594.5234375 ... 102594.5234375
  102594.5234375 102594.5234375]
 [102586.5234375 102586.5234375 102586.5234375 ... 102586.5234375
  102586.5234375 102586.5234375]
 [102574.5234375 102574.5234375 102574.5234375 ... 102574.5234375
  102574.5234375 102574.5234375]]
Writing has been done. Trying to read.
ECCODES ERROR   :  Data section size mismatch: offset before data=175, offset after data=379002 (num values=1038240, bits per value=12)
Traceback (most recent call last):
  File "/home/chabotv/tmp/test_epy.py", line 21, in <module>
    MSL = resAearo.readfield({
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 2559, in readfield
    matchingfields = self.readfields(handgrip,
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/resources/FileResource.py", line 51, in nowopen
    return mtd(self, *args, **kwargs)
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 2636, in readfields
    matchingfields.append(msg.as_field(getdata=getdata,
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 1571, in as_field
    field.setdata(self._read_data(field.geometry))
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 1903, in _read_data
    data1d = self['values']
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 502, in __getitem__
    self._readattribute(item)
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 648, in _readattribute
    raise e
  File "/home/chabotv/.local/lib/python3.10/site-packages/EPyGrAM/epygram/formats/GRIB.py", line 645, in _readattribute
    attr = lowlevelgrib.get_values(self._gid)
  File "/usr/lib/python3/dist-packages/gribapi/gribapi.py", line 2005, in grib_get_values
    return grib_get_double_array(gribid, "values")
  File "/usr/lib/python3/dist-packages/gribapi/gribapi.py", line 1173, in grib_get_double_array
    GRIB_CHECK(err)
  File "/usr/lib/python3/dist-packages/gribapi/gribapi.py", line 228, in GRIB_CHECK
    errors.raise_grib_error(errid)
  File "/usr/lib/python3/dist-packages/gribapi/errors.py", line 382, in raise_grib_error
    raise ERROR_MAP[errid](errid)
gribapi.errors.DecodingError: Decoding invalid
Grib dump before modification ***** FILE: ARPEGE_2023100612_MER_GLOB025.grib #============== MESSAGE 1 ( length=1557539 ) ============== GRIB { # Meteorological products (grib2/tables/15/0.0.table) discipline = 0; editionNumber = 2; # French Weather Service - Toulouse (common/c-11.table) centre = 85; subCentre = 0; # Start of forecast (grib2/tables/15/1.2.table) significanceOfReferenceTime = 1; dataDate = 20231006; dataTime = 1200; # Operational products (grib2/tables/15/1.3.table) productionStatusOfProcessedData = 0; # Forecast products (grib2/tables/15/1.4.table) typeOfProcessedData = 1; numberOfDataPoints = 1038240; # There is no appended list (grib2/tables/15/3.11.table) interpretationOfNumberOfPoints = 0; # Latitude/longitude (Also called equidistant cylindrical, or Plate Carree) (grib2/tables/15/3.1.table) gridDefinitionTemplateNumber = 0; # Earth assumed spherical with radius of 6 371 229.0 m (grib2/tables/15/3.2.table) shapeOfTheEarth = 6; Ni = 1440; Nj = 721; iScansNegatively = 0; jScansPositively = 0; jPointsAreConsecutive = 0; alternativeRowScanning = 0; latitudeOfFirstGridPointInDegrees = 90; longitudeOfFirstGridPointInDegrees = 0; latitudeOfLastGridPointInDegrees = -90; longitudeOfLastGridPointInDegrees = 359.75; iDirectionIncrementInDegrees = 0.25; jDirectionIncrementInDegrees = 0.25; gridType = regular_ll; NV = 0; # Analysis or forecast at a horizontal level or in a horizontal layer at a point in time (grib2/tables/15/4.0.table) productDefinitionTemplateNumber = 0; # Mass (grib2/tables/15/4.1.0.table) parameterCategory = 3; # Pressure reduced to MSL (Pa) (grib2/tables/15/4.2.0.3.table) parameterNumber = 1; #-READ ONLY- parameterUnits = Pa; #-READ ONLY- parameterName = Pressure reduced to MSL ; # Forecast (grib2/tables/15/4.3.table) typeOfGeneratingProcess = 2; generatingProcessIdentifier = 211; # Hour (grib2/tables/15/4.4.table) indicatorOfUnitOfTimeRange = 1; # Hour (stepUnits.table) stepUnits = 1; forecastTime = 0; stepRange = 0; # Mean sea level (grib2/tables/15/4.5.table) typeOfFirstFixedSurface = 101; #-READ ONLY- unitsOfFirstFixedSurface = unknown; #-READ ONLY- nameOfFirstFixedSurface = Mean sea level ; scaleFactorOfFirstFixedSurface = 0; scaledValueOfFirstFixedSurface = 0; # Missing (grib2/tables/15/4.5.table) typeOfSecondFixedSurface = 255; #-READ ONLY- unitsOfSecondFixedSurface = unknown; #-READ ONLY- nameOfSecondFixedSurface = Missing ; scaleFactorOfSecondFixedSurface = MISSING; scaledValueOfSecondFixedSurface = MISSING; level = 0; shortNameLegacyECMF = unknown; shortNameECMF = prmsl; shortName = prmsl; nameLegacyECMF = unknown; nameECMF = Pressure reduced to MSL; name = Pressure reduced to MSL; cfNameLegacyECMF = unknown; cfNameECMF = unknown; cfName = unknown; cfVarNameLegacyECMF = unknown; cfVarNameECMF = prmsl; cfVarName = prmsl; #-READ ONLY- modelName = unknown; numberOfValues = 1038240; packingType = grid_simple; # A bit map does not apply to this product (grib2/tables/15/6.0.table) bitMapIndicator = 255; bitmapPresent = 0; values(1038240) = { 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575, 102575 ... 1038140 more values } #-READ ONLY- maximum = 104339; #-READ ONLY- minimum = 94894.5; #-READ ONLY- average = 100825; #-READ ONLY- numberOfMissing = 0; #-READ ONLY- standardDeviation = 1533.86; #-READ ONLY- skewness = -1.09955; #-READ ONLY- kurtosis = 0.92373; #-READ ONLY- isConstant = 0; #-READ ONLY- getNumberOfValues = 1038240; }
AlexandreMary commented 9 months ago

There seem to be an inconsistency between the packingType of the sample used (by default) to generate the GRIB message and the packing type requested ('grid_second_order', by default). A workaround (before I fix that) is to specify the packing type at write time: resAearo.writefield(MSL, other_GRIB_options={'packingType':'grid_simple'})