grib.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Mar 15, 02:47

grib.c
View Options

	/*
	* Copyright (C) iMariner. All Rights Reserved.
	* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER
	*
	* This program is free software; you can redistribute it and/or
	* modify it under the terms of the GNU General Public License version
	* 2 only, as published by the Free Software Foundation.
	*
	* This program is distributed in the hope that it will be useful, but
	* WITHOUT ANY WARRANTY; without even the implied warranty of
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	* General Public License version 2 for more details (a copy is
	* included at /legal/license.txt).
	*
	* You should have received a copy of the GNU General Public License
	* version 2 along with this work; if not, write to the Free Software
	* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
	* 02110-1301 USA
	*
	*/

	/**
	Extracts from ZyGrib / GribV1Record.cpp
	*/

	/**
	World Meteorogical Organization WMO-No.306
	Manual on Codes
	International Codes
	VOLUME I.2
	PART B - Binary Codes
	http://library.wmo.int/opac/index.php?lvl=notice_display&id=10684
	https://www.wmo.int/pages/prog/www/WMOCodes/WMO306_vI2/Publications/2011editionUP2014/WMO306_vI2_2011UP2014.pdf
	*/

	#include <stdio.h>
	#include <stdlib.h>
	#include <stdint.h>
	#include <stddef.h>
	#include <stdbool.h>
	#include <sys/types.h>
	#include <unistd.h>
	#include <string.h>
	#include <math.h>
	#include <fcntl.h>
	#include <time.h>


	#define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
	#define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))

	// Table Code 2
	#define DATATYPE_PRESSURE 2
	#define DATATYPE_TEMP 11 // K
	#define DATATYPE_UWIND 33 // m/s
	#define DATATYPE_VWIND 34 // m/s
	#define DATATYPE_PRECIP 61 // kg/m2
	#define DATATYPE_WAVES 100 // m
	#define DATATYPE_APCP 157

	// Table Code 3
	#define LEVELTYPE_NO 0
	#define LEVELTYPE_GROUND 1
	#define LEVELTYPE_MSL 102 // Mean Sea Level
	#define LEVELTYPE_HEIGHT 105


	bool Verbose = false;

	unsigned char buf[1024*100];

	// IS : Indicator Section
	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize;
	unsigned long totalLength;
	unsigned char editionNumber; // GRIB1 or GRIB2
	char magic[5];
	} T_section0;

	// PDS : Product Description Section
	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize; // Length of section
	unsigned char parameterTableVersion; // GRIB tables Version No
	unsigned char idCenter;
	unsigned char idModel;
	unsigned char idGrid;
	unsigned char flags;
	unsigned char dataType; // PDS9
	unsigned char levelType; // PDS10
	unsigned int levelValue; // PDS11 & PDS12
	unsigned int year;
	unsigned char month;
	unsigned char day;
	unsigned char hour;
	unsigned char minute;
	unsigned char forcastTimeUnit;
	unsigned char periodP1;
	unsigned char periodP2;
	unsigned char timeRange;
	int D;
	double Dpow10;
	} T_section1;

	// GDS : Grid Description Section
	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize;
	unsigned char NV; // number of vertical coordinate parameters
	unsigned char PV; // PV - location (octet number) of the list of numbers
	// of points in each row (if no vertical coordinate
	// parameters are present), if present;
	// or PL – location (octet number) of the list of numbers
	// of points in each row (if no vertical coordinate
	// parameters are present), if present;
	// or 255 (all bits set to 1) if neither are present
	unsigned char gridType; // Data representation type (see Code table 6)
	unsigned int Ni;
	unsigned int Nj;
	double La1;
	double Lo1;
	unsigned char resolFlags;
	double La2;
	double Lo2;
	double Di;
	double Dj;
	unsigned char scanFlags;
	double lonMin;
	double lonMax;
	double latMin;
	double latMax;

	bool hasDiDj;
	bool isEarthSpheric;
	bool isUeastVnorth;
	bool isScanIpositive;
	bool isScanJpositive;
	bool isAdjacentI;
	} T_section2;

	// BMS : BitMap Section
	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize;
	} T_section3;

	// BDS : Binary Data Section
	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize;
	unsigned char flags;
	int scaleFactorE;
	double refValue;
	unsigned char nbBitsInPack;
	double scaleFactorEpow2;
	unsigned char unusedBitsEndBDS;

	bool isGridData;
	bool isSimplePacking;
	bool isFloatValues;
	double data[100*1024];
	} T_section4;

	// ES : End Section
	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize;
	char magic[5];
	} T_section5;

	typedef struct {
	T_section0 s0;
	T_section1 s1;
	T_section2 s2;
	T_section3 s3;
	T_section4 s4;
	T_section5 s5;
	} T_gribRecord;

	//T_gribRecord rec;

	unsigned int readInt2(unsigned char *buf) {
	return buf[0]<<8 \| buf[1];
	}

	unsigned int readInt3(unsigned char *buf) {
	return buf[0]<<16 \| buf[1]<<8 \| buf[2];
	}

	unsigned int readInt4(unsigned char *buf) {
	return buf[0]<<24 \| buf[1]<<16 \| buf[2]<<8 \| buf[3];
	}

	unsigned long readLong8(unsigned char *buf) {
	int i;
	unsigned long l = 0;

	for (i=0; i<8; i++)
	l = (l<<8) + buf[i];
	return l;
	}

	int readSignedInt2(unsigned char *buf) {
	int sign = buf[0]&0x80 ? -1 : 1;
	return ((buf[0]&0x7f)<<8 \| buf[1]) * sign;
	}

	int readSignedInt3(unsigned char *buf) {
	int sign = buf[0]&0x80 ? -1 : 1;
	return ((buf[0]&0x7f)<<16 \| buf[1]<<8 \| buf[2])*sign;
	}

	double readFloat4(unsigned char *buf) {
	double val;
	int sign = buf[0]&0x80 ? -1 : 1;
	int A = buf[0]&0x7F;
	int B = readInt3(buf+1);

	val = pow(2.,-24)Bpow(16.,A-64);
	return val * sign;
	}


	/* wesley ebisuzaki v0.2
	*
	* takes 4 byte character string (single precision ieee big-endian)
	* and returns a float
	*
	* doesn't handle NAN, infinity and any other funny stuff in ieee
	*
	* ansi C
	*/

	float ieee2flt(unsigned char *ieee) {
	double fmant;
	int exp;

	if ((ieee[0] & 127) == 0 && ieee[1] == 0 && ieee[2] == 0 && ieee[3] == 0)
	return (float) 0.0;

	exp = ((ieee[0] & 127) << 1) + (ieee[1] >> 7);
	fmant = (double) ((int) ieee[3] + (int) (ieee[2] << 8) +
	(int) ((ieee[1] \| 128) << 16));
	if (ieee[0] & 128) fmant = -fmant;
	return (float) (ldexp(fmant, (int) (exp - 128 - 22)));
	}

	unsigned int readPackedBits(unsigned char *buf, unsigned int first, unsigned int nbBits) {
	unsigned int oct = first / 8;
	unsigned int bit = first % 8;
	//printf ("readPackedBits(%d,%d)\n", first, nbBits);
	unsigned int val = (buf[oct]<<24) + (buf[oct+1]<<16) + (buf[oct+2]<<8) + (buf[oct+3]);
	val = val << bit;
	val = val >> (32-nbBits);
	return val;
	}

	// PDS : Product Definition Section
	int GRIB1_readSection1(int fd, T_gribRecord *rec) {
	off_t offset = rec->s0.offsetInFile+rec->s0.sectionSize;
	if (lseek (fd, offset, SEEK_SET) != offset)
	return -1;
	rec->s1.offsetInFile = offset;
	//printf ("s1.offsetInFile = %ld\n", rec->s1.offsetInFile);

	if (read (fd, buf, 3) != 3)
	return -2;
	rec->s1.sectionSize = readInt3(buf);
	//printf ("s1.sectionSize = %ld\n", rec->s1.sectionSize);

	if (read (fd, buf, rec->s1.sectionSize-3) != rec->s1.sectionSize-3)
	return -3;

	rec->s1.parameterTableVersion = buf[0];
	rec->s1.idCenter = buf[1];
	rec->s1.idModel = buf[2];
	rec->s1.idGrid = buf[3];
	rec->s1.flags = buf[4];
	rec->s1.dataType = buf[5];
	rec->s1.levelType = buf[6];
	rec->s1.levelValue = readInt2(buf+7);
	rec->s1.year = (buf[21]-1)*100+buf[9];
	rec->s1.month = buf[10];
	rec->s1.day = buf[11];
	rec->s1.hour = buf[12];
	rec->s1.minute = buf[13];
	rec->s1.forcastTimeUnit = buf[14];
	rec->s1.timeRange= buf[17];
	switch (rec->s1.timeRange) {
	case 10 :
	rec->s1.periodP1 = (buf[15] << 8) + buf[16];
	rec->s1.periodP2 = 0;
	break;
	default :
	rec->s1.periodP1 = buf[15];
	rec->s1.periodP2 = buf[16];
	break;
	}

	int decim = (buf[23]&0x7f)<<8 \| buf[24];
	if (buf[23] & 0x80)
	decim = -decim;
	rec->s1.D = decim;
	rec->s1.Dpow10 = pow(10.0, decim);
	return 0;
	}

	void GRIB1_dumpSection1(T_gribRecord *rec) {
	if (!Verbose) return;
	printf ("PDS :\n");
	printf ("s1.offsetInFile = %ld\n", rec->s1.offsetInFile);
	printf ("s1.sectionSize = %ld\n", rec->s1.sectionSize);
	printf ("s1.parameterTableVersion = %d : GRIB tables Version No. (currently 3 for international exchange)\n", rec->s1.parameterTableVersion);
	printf ("s1.idCenter = %d : Identification of originating/generating centre (see Code table 0 = Common Code table C–1 in Part C/c.)\n", rec->s1.idCenter);
	printf ("s1.idModel = %d\n", rec->s1.idModel);
	printf ("s1.idGrid = %d\n", rec->s1.idGrid);
	printf ("s1.flags = 0x%02x\n", rec->s1.flags);
	printf ("s1.dataType = %d : Indicator of parameter (see Code table 2) 33=u-wind, 34=v-wind, 2=pressure\n", rec->s1.dataType);
	printf ("s1.levelType = %d : Indicator of type of level (see Code table 3)\n", rec->s1.levelType);
	printf ("s1.levelValue = %d : Height, pressure, etc. of levels (see Code table 3)\n", rec->s1.levelValue);
	printf ("s1.year = %d : Reference time of data - date and time of start of averaging or accumulation period\n", rec->s1.year);
	printf ("s1.month = %d\n", rec->s1.month);
	printf ("s1.day = %d\n", rec->s1.day);
	printf ("s1.hour = %d\n", rec->s1.hour);
	printf ("s1.minute = %d\n", rec->s1.minute);
	printf ("s1.forcastTimeUnit = %d : Indicator of unit of time range (see Code table 4) 0=min 1=hour 2=day ...\n", rec->s1.forcastTimeUnit);
	printf ("s1.periodP1 = %d : Period of time (number of time units) (0 for analyses or initialized analyses). Units of time given by octet 18\n", rec->s1.periodP1);
	printf ("s1.periodP2 = %d : Period of time (number of time units); or Time interval between successive analyses, initialized analyses or forecasts, undergoing averaging or accumulation. Units of time given by octet 18\n", rec->s1.periodP2);
	printf ("s1.timeRange = %d : Time range indicator (see Code table 5)\n", rec->s1.timeRange);
	printf ("s1.D = %d : Units decimal scale factor (D)\n", rec->s1.D);
	printf ("s1.Dpow10 = %f\n", rec->s1.Dpow10);
	printf ("\n");
	}

	int GRIB1_readSection2(int fd, T_gribRecord *rec) {
	// There is no GDS
	if ((rec->s1.flags & 0x80) == 0) {
	rec->s2.offsetInFile = 0;
	rec->s2.sectionSize = 0;
	return 0;
	}
	off_t offset = rec->s1.offsetInFile+rec->s1.sectionSize;
	if (lseek (fd, offset, SEEK_SET) != offset) {
	printf ("ERROR LSEEK %ld\n", offset);
	return -1;
	}
	rec->s2.offsetInFile = offset;

	if (read (fd, buf, 3) != 3)
	return -2;
	rec->s2.sectionSize = readInt3(buf);

	if (read (fd, buf, rec->s2.sectionSize-3) != rec->s2.sectionSize-3)
	return -3;

	rec->s2.NV = buf[0];
	rec->s2.PV = buf[1];
	rec->s2.gridType = buf[2];

	// It is not a Latitude/Longitude Grid
	if (rec->s2.gridType != 0) {
	printf ("FATAL ERROR : It is not a Latitude/Longitude Grid\n");
	return -1;
	}

	rec->s2.Ni = readInt2(buf+3);
	rec->s2.Nj = readInt2(buf+5);
	rec->s2.La1 = readSignedInt3(buf+7)/1000.0;
	rec->s2.Lo1 = readSignedInt3(buf+10)/1000.0;
	rec->s2.resolFlags = buf[13];
	rec->s2.La2 = readSignedInt3(buf+14)/1000.0;
	rec->s2.Lo2 = readSignedInt3(buf+17)/1000.0;

	if (rec->s2.Lo1>=0 && rec->s2.Lo1<=180 && rec->s2.Lo2<0) {
	rec->s2.Lo2 += 360.0; // cross the 180 deg meridien,beetwen alaska and russia
	}

	rec->s2.Di = readSignedInt2(buf+20)/1000.0;
	rec->s2.Dj = readSignedInt2(buf+22)/1000.0;

	while ( rec->s2.Lo1> rec->s2.Lo2 && rec->s2.Di >0) { // horizontal size > 360 °
	rec->s2.Lo1 -= 360.0;
	}

	rec->s2.hasDiDj = (rec->s2.resolFlags&0x80) !=0;
	rec->s2.isEarthSpheric = (rec->s2.resolFlags&0x40) ==0;
	rec->s2.isUeastVnorth = (rec->s2.resolFlags&0x08) ==0;

	rec->s2.scanFlags = buf[24];

	rec->s2.isScanIpositive = (rec->s2.scanFlags&0x80) ==0;
	rec->s2.isScanJpositive = (rec->s2.scanFlags&0x40) !=0;
	rec->s2.isAdjacentI = (rec->s2.scanFlags&0x20) ==0;

	if (rec->s2.Lo2 > rec->s2.Lo1) {
	rec->s2.lonMin = rec->s2.Lo1;
	rec->s2.lonMax = rec->s2.Lo2;
	}
	else {
	rec->s2.lonMin = rec->s2.Lo2;
	rec->s2.lonMax = rec->s2.Lo1;
	}
	if (rec->s2.La2 > rec->s2.La1) {
	rec->s2.latMin = rec->s2.La1;
	rec->s2.latMax = rec->s2.La2;
	}
	else {
	rec->s2.latMin = rec->s2.La2;
	rec->s2.latMax = rec->s2.La1;
	}
	if (rec->s2.Ni<=1 \|\| rec->s2.Nj<=1) {
	return -1;
	}
	else {
	rec->s2.Di = (rec->s2.Lo2-rec->s2.Lo1) / (rec->s2.Ni-1);
	rec->s2.Dj = (rec->s2.La2-rec->s2.La1) / (rec->s2.Nj-1);
	}

	return 0;
	}

	void GRIB1_dumpSection2(T_gribRecord *rec) {
	if (!Verbose) return;
	printf ("GDS :\n");
	printf ("s2.offsetInFile = %ld\n", rec->s2.offsetInFile);
	printf ("s2.sectionSize = %ld\n", rec->s2.sectionSize);
	printf ("s2.NV = %d : number of vertical coordinate parameters\n", rec->s2.NV);
	printf ("s2.PV = %d\n", rec->s2.PV);
	printf ("s2.gridType = %d : Data representation type (see Code table 6). 0 in our case (Latitude/longitude).\n", rec->s2.gridType);
	printf ("s2.Ni = %d (number of points along a parallel)\n", rec->s2.Ni);
	printf ("s2.Nj = %d (number of points along a meridian)\n", rec->s2.Nj);
	printf ("s2.La1 = %f (latitude of first grid point)\n", rec->s2.La1);
	printf ("s2.Lo1 = %f (longitude of first grid point)\n", rec->s2.Lo1);
	printf ("s2.resolFlags = %x : Resolution and component flags (see Code table 7)\n", rec->s2.resolFlags);
	printf (" %s %s %s\n",
	rec->s2.hasDiDj ? "hasDiDj":"",
	rec->s2.isEarthSpheric ? "isEarthSpheric":"",
	rec->s2.isUeastVnorth ? "isUeastVnorth":"");
	printf ("s2.La2 = %f (latitude of last grid point)\n", rec->s2.La2);
	printf ("s2.Lo2 = %f (longitude of last grid point)\n", rec->s2.Lo2);
	printf ("s2.Di = %f (i direction increment)\n", rec->s2.Di);
	printf ("s2.Dj = %f (j direction increment)\n", rec->s2.Dj);
	printf ("s2.scanFlags = %x : Scanning mode (flags – see Flag/Code table 8)\n", rec->s2.scanFlags);
	if (rec->s2.scanFlags & 0x80)
	printf (" Points scan in -i direction (east to west)\n");
	else
	printf (" Points scan in +i direction (west to east)\n");
	if (rec->s2.scanFlags & 0x40)
	printf (" Points scan in +j direction (south to north)\n");
	else
	printf (" Points scan in -j direction (north to south)\n");
	if (rec->s2.scanFlags & 0x20)
	printf (" Adjacent points in j direction are consecutive\n");
	else
	printf (" Adjacent points in i direction are consecutive\n");
	printf ("\n");
	}

	int GRIB1_readSection3(int fd, T_gribRecord *rec) {
	off_t offset;
	// There is no BMS
	if ((rec->s1.flags & 0x40) == 0) {
	rec->s3.offsetInFile = 0;
	rec->s3.sectionSize = 0;
	return -1;
	}
	if (rec->s1.flags & 0x80)
	offset = rec->s2.offsetInFile+rec->s2.sectionSize;
	else
	offset = rec->s1.offsetInFile+rec->s1.sectionSize;
	if (lseek (fd, offset, SEEK_SET) != offset)
	return -1;
	rec->s3.offsetInFile = offset;

	if (read (fd, buf, 3) != 3)
	return -2;
	rec->s3.sectionSize = readInt3(buf);

	if (read (fd, buf, rec->s3.sectionSize-3) != rec->s3.sectionSize-3)
	return -3;

	// TBD

	return 0;
	}

	void GRIB1_dumpSection3(T_gribRecord *rec) {
	if (!Verbose) return;
	printf ("BMS :\n");
	printf ("s3.offsetInFile = %ld\n", rec->s3.offsetInFile);
	printf ("s3.sectionSize = %ld\n", rec->s3.sectionSize);
	printf ("\n");
	}

	void GRIB1_dumpSection4(T_gribRecord *rec) {
	if (!Verbose) return;
	printf ("BDS :\n");
	printf ("s4.offsetInFile = %ld\n", rec->s4.offsetInFile);
	printf ("s4.sectionSize = %ld\n", rec->s4.sectionSize);

	printf ("s4.flags = %x\n", rec->s4.flags);
	printf ("s4.unusedBitsEndBDS = %d\n", rec->s4.unusedBitsEndBDS);
	printf ("s4.scaleFactorE = %d\n", rec->s4.scaleFactorE);
	printf ("s4.refValue = %f : minimum of packed values\n", rec->s4.refValue);
	printf ("s4.nbBitsInPack = %d : Number of bits containing each packed value\n", rec->s4.nbBitsInPack);
	printf ("s4.scaleFactorEpow2 = %f\n", rec->s4.scaleFactorEpow2);
	printf ("\n");
	}

	int GRIB1_readSection4(int fd, T_gribRecord *rec) {
	off_t offset;
	if (rec->s1.flags & 0x40)
	offset = rec->s3.offsetInFile+rec->s3.sectionSize;
	else if (rec->s1.flags & 0x80)
	offset = rec->s2.offsetInFile+rec->s2.sectionSize;
	else
	offset = rec->s1.offsetInFile+rec->s1.sectionSize;

	rec->s4.offsetInFile = offset;

	if (lseek (fd, offset, SEEK_SET) != offset)
	return -1;
	if (read (fd, buf, 3) != 3)
	return -2;
	rec->s4.sectionSize = readInt3(buf);

	int n = read (fd, buf, rec->s4.sectionSize-3);
	if (n != rec->s4.sectionSize-3)
	return -3;

	rec->s4.flags = buf[0] & 0xF0;
	rec->s4.unusedBitsEndBDS = buf[0] & 0x0F;
	rec->s4.scaleFactorE = readSignedInt2(buf+1);
	rec->s4.refValue = readFloat4(buf+3);
	rec->s4.nbBitsInPack = buf[7];
	rec->s4.scaleFactorEpow2 = pow(2.,rec->s4.scaleFactorE);
	rec->s4.isGridData = (rec->s4.flags&0x80) ==0;
	rec->s4.isSimplePacking = (rec->s4.flags&0x80) ==0;
	rec->s4.isFloatValues = (rec->s4.flags&0x80) ==0;

	GRIB1_dumpSection4(rec);

	// Read data in the order given by isAdjacentI
	int i, j, x;
	int ind;
	unsigned int startbit = 0;

	if (rec->s2.isAdjacentI) {
	//double La = MIN(rec->s2.La1, rec->s2.La2);
	double La = rec->s2.La1;
	for (j=0; j<rec->s2.Nj; j++) {
	//double Lo = MIN(rec->s2.Lo1, rec->s2.Lo2);
	double Lo = rec->s2.Lo1;
	for (i=0; i<rec->s2.Ni; i++) {
	if (!rec->s2.hasDiDj && !rec->s2.isScanJpositive) {
	//if (!rec->s2.isScanJpositive) {
	ind = (rec->s2.Nj-1 -j)*rec->s2.Ni+i;
	}
	else {
	ind = j*rec->s2.Ni+i;
	}
	// if (rec->s2.hasValue(i,j)) {
	x = readPackedBits(buf+8, startbit, rec->s4.nbBitsInPack);
	rec->s4.data[ind] = (rec->s4.refValue + x*rec->s4.scaleFactorEpow2)/rec->s1.Dpow10;
	if (Verbose)
	printf ("La:%f Lo:%f ind:%d data[%d,%d] = %d => %f\n", La, Lo, ind, i, j, x, rec->s4.data[ind]);
	startbit += rec->s4.nbBitsInPack;
	/*
	}
	else {
	rec->s4.data[ind] = GRIB_NOTDEF;
	}
	*/
	Lo = Lo + rec->s2.Di;
	}
	La = La + rec->s2.Dj;
	}
	}
	else {
	for (i=0; i<rec->s2.Ni; i++) {
	for (j=0; j<rec->s2.Nj; j++) {
	if (!rec->s2.hasDiDj && !rec->s2.isScanJpositive) {
	ind = (rec->s2.Nj-1 -j)*rec->s2.Ni+i;
	}
	else {
	ind = j*rec->s2.Ni+i;
	}
	// if (rec->s2.hasValue(i,j)) {
	x = readPackedBits(buf+8, startbit, rec->s4.nbBitsInPack);
	startbit += rec->s4.nbBitsInPack;
	rec->s4.data[ind] = (rec->s4.refValue + x*rec->s4.scaleFactorEpow2)/rec->s1.Dpow10;
	//printf ("rec->s4.data[%d,%d] = %f\n", i, j, rec->s4.data[ind]);
	/*
	}
	else {
	rec->s4.data[ind] = GRIB_NOTDEF;
	}
	*/
	}
	}
	}

	return 0;
	}

	double getDataAtIndex(T_gribRecord *rec, int i, int j) {
	int ind = j*rec->s2.Ni+i;
	return rec->s4.data[ind];
	}

	double getDataAtValues(T_gribRecord *rec, double La, double Lo) {
	int i = (Lo - rec->s2.Lo1) / rec->s2.Di;
	int j = (La - rec->s2.La1) / rec->s2.Dj;
	if (Verbose)
	printf ("getDataAtIndex (%f,%f) (i:%d,j:%d)\n", La, Lo, i, j);
	return getDataAtIndex (rec, i, j);
	}

	bool inRectangle(T_gribRecord *rec, double La, double Lo) {
	//printf ("XXX %f,%f %f %f %f %f\n", La, Lo, MIN(rec->s2.La1, rec->s2.La2), MAX(rec->s2.La1, rec->s2.La2), MIN(rec->s2.Lo1, rec->s2.Lo2), MAX(rec->s2.Lo1, rec->s2.Lo2));
	if (La < MIN(rec->s2.La1, rec->s2.La2))
	return false;
	if (La > MAX(rec->s2.La1, rec->s2.La2))
	return false;
	if (Lo < MIN(rec->s2.Lo1, rec->s2.Lo2))
	return false;
	if (Lo > MAX(rec->s2.Lo1, rec->s2.Lo2))
	return false;
	return true;
	}

	int GRIB1_readSection5(int fd, T_gribRecord *rec) {
	off_t offset = rec->s4.offsetInFile+rec->s4.sectionSize;
	if (lseek (fd, offset, SEEK_SET) != offset)
	return -1;
	rec->s5.offsetInFile = offset;
	rec->s5.sectionSize = 4;

	if (read (fd, buf, rec->s5.sectionSize) != rec->s5.sectionSize)
	return -1;
	memcpy (rec->s5.magic, buf, 4);
	rec->s5.magic[4] = 0;
	if (strcmp(rec->s5.magic, "7777"))
	return -2;
	return 0;
	}

	void GRIB1_dumpSection5( T_gribRecord *rec) {
	if (!Verbose) return;
	printf ("ES :\n");
	printf ("s5.offsetInFile = %ld\n", rec->s5.offsetInFile);
	printf ("s5.sectionSize = %ld\n", rec->s5.sectionSize);
	printf ("s5.magic = %s\n", rec->s5.magic);
	printf ("\n");
	}

	int GRIB1_seekToNextRecord(int fd, T_gribRecord *rec) {
	off_t offset = rec->s0.offsetInFile+rec->s0.totalLength;
	off_t res;
	res = lseek (fd, offset, SEEK_SET);
	//printf ("seekToNextRecord nextRecord %ld. res=%ld\n", offset, res);
	return res;
	}

	int GRIB1(int fd, T_gribRecord *rec) {
	GRIB1_readSection1 (fd, rec);
	GRIB1_dumpSection1 (rec);

	GRIB1_readSection2 (fd, rec);
	GRIB1_dumpSection2 (rec);

	GRIB1_readSection3 (fd, rec);
	GRIB1_dumpSection3 (rec);

	GRIB1_readSection4 (fd, rec);
	//GRIB1_dumpSection4 (rec);

	GRIB1_readSection5 (fd, rec);
	GRIB1_dumpSection5 (rec);
	return 0;
	}





	//====================================== GRIB2 =============================

	typedef struct {
	off_t offsetInFile;
	unsigned long sectionSize;
	unsigned long totalLength;
	unsigned char editionNumber;
	char magic[5];
	unsigned int discipline;
	} T_grib2sect0;

	typedef struct {
	unsigned int center;
	unsigned int subCenter;
	unsigned int year;
	unsigned char month;
	unsigned char day;
	unsigned char hour;
	unsigned char minute;
	} T_grib2sect1;

	typedef struct {
	int Ni, Nj;
	double La1, Lo1, La2, Lo2;
	double Di, Dj;
	} T_grib2sect3;

	typedef struct {
	int E;
	int D;
	double R;
	double Epow2;
	double Dpow10;
	int nbBits;
	} T_grib2sect5;

	typedef struct {
	double data[100*1024];
	} T_grib2sect7;

	typedef struct {
	T_grib2sect0 s0;
	T_grib2sect1 s1;
	T_grib2sect3 s3;
	T_grib2sect5 s5;
	T_grib2sect7 s7;
	} T_grib2Record;

	//T_grib2Record rec2;

	int GRIB2_readSection1 (int fd, T_grib2Record *rec2, int secLen) {
	if (read (fd, buf, secLen) != secLen)
	return -1;

	rec2->s1.center = readInt2(buf);
	rec2->s1.subCenter = readInt2(buf+2);

	printf ("==== Identification Section====\n");
	printf ("rec2->s1.center = %d\n", rec2->s1.center);
	printf ("rec2->s1.subCenter = %d\n", rec2->s1.subCenter);
	printf ("rec2->s1.masterTables = %d\n", buf[4]);
	printf ("Significance of ref time = %d (0:analysis 1:start of forecast)\n", buf[6]);
	printf ("Year = %d\n", readInt2(buf+7));
	printf ("Month = %d\n", buf[9]);
	printf ("Day = %d\n", buf[10]);
	printf ("Hour = %d\n", buf[11]);
	printf ("Minute = %d\n", buf[12]);
	printf ("Second = %d\n", buf[13]);
	printf ("Production Status = %d\n", buf[14]);
	printf ("Type of processed data = %d (0:analysis 1:forecast)\n", buf[15]);
	return 0;
	}

	int GRIB2_readSection3 (int fd, T_grib2Record *rec2, int secLen) {
	if (read (fd, buf, secLen) != secLen)
	return -1;

	printf ("==== Grid Definition Section ====\n");
	printf ("Number of data points = %d\n", readInt4(buf+1));
	printf ("Grid template = %d\n", readInt2(buf+7));
	switch (readInt2(buf+7)) {
	case 0: // Lat/Lon
	printf ("Latitude/Longitude\n");
	rec2->s3.Ni = readInt4(buf+25);
	rec2->s3.Nj = readInt4(buf+29);
	rec2->s3.La1 = readInt4(buf+41)/1000000.0;
	rec2->s3.Lo1 = readInt4(buf+45)/1000000.0;
	rec2->s3.La2 = readInt4(buf+50)/1000000.0;
	rec2->s3.Lo2 = readInt4(buf+54)/1000000.0;
	rec2->s3.Di = readInt4(buf+58)/1000000.0;
	rec2->s3.Dj = readInt4(buf+62)/1000000.0;

	printf ("Ni = %d\n", rec2->s3.Ni);
	printf ("Nj = %d\n", rec2->s3.Nj);
	printf ("La1 = %f\n", rec2->s3.La1);
	printf ("Lo1 = %f\n", rec2->s3.Lo1);
	printf ("La2 = %f\n", rec2->s3.La2);
	printf ("Lo2 = %f\n", rec2->s3.Lo2);
	printf ("Di = %f\n", rec2->s3.Di);
	printf ("Dj = %f\n", rec2->s3.Di);
	break;
	}
	return 0;
	}

	int GRIB2_readSection4 (int fd, T_grib2Record *rec2, int secLen) {
	if (read (fd, buf, secLen) != secLen)
	return -1;

	printf ("==== Product Definition Section ====\n");
	printf ("Product Definition Template = %d\n", readInt2(buf+2));
	switch (readInt2(buf+2)) {
	case 0: // Analysis or forecast at a horizontal level
	printf ("Analysis or forecast at a horizontal level\n");
	printf ("Parameter category = %d\n", buf[4]);
	printf ("Parameter number = %d\n", buf[5]);
	switch (rec2->s0.discipline) {
	case 0 : // Meteorological Products
	// http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-1.shtml#0
	switch (buf[4]) {
	case 2 : // Momentum
	if (buf[5] == 0)
	printf ("=> WDIR (°)\n");
	else if (buf[5] == 1)
	printf ("=> WIND (m/s)\n");
	else if (buf[5] == 2)
	printf ("=> UGRD (m/s)\n");
	else if (buf[5] == 3)
	printf ("=> VGRD (m/s)\n");
	break;
	}
	break;
	case 10 : // Oceanographic Products
	switch (buf[4]) {
	case 0 : // Waves
	// http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-2-10-0.shtml
	if (buf[5] == 3)
	printf ("=> HTSGW (m)\n");
	else if (buf[5] == 4)
	printf ("=> WVDIR (Degree True)\n");
	else if (buf[5] == 5)
	printf ("=> WVHGT (m)\n");
	else if (buf[5] == 6)
	printf ("=> WVPER (s)\n");
	else if (buf[5] == 7)
	printf ("=> SWDIR (Degree True)\n");
	else if (buf[5] == 8)
	printf ("=> SWELL (m)\n");
	else if (buf[5] == 9)
	printf ("=> WSPER (s)\n");
	else if (buf[5] == 10)
	printf ("=> DIRPW (Degree True)\n");
	else if (buf[5] == 11)
	printf ("=> PERPW (s)\n");
	else if (buf[5] == 12)
	printf ("=> DIRSW (Degree True)\n");
	else if (buf[5] == 13)
	printf ("=> PERSW (s)\n");
	else if (buf[5] == 14)
	printf ("=> WWSDIR (Degree True)\n");
	else if (buf[5] == 15)
	printf ("=> MWSPER (s)\n");
	break;
	case 1 : // Currents
	if (buf[5] == 0)
	printf ("=> DIR C (Degree True)\n");
	if (buf[5] == 1)
	printf ("=> SP C (m/s)\n");
	if (buf[5] == 2)
	printf ("=> UOGRD (m/s)\n");
	if (buf[5] == 3)
	printf ("=> VOGRD (m/s)\n");
	break;
	}
	break;
	}
	printf ("Hours after ref time = %d\n", readInt2(buf+9));
	printf ("Minutes after ref time = %d\n", buf[11]);
	printf ("Indicator of units of time range : %d (1:hour 2:day)\n", buf[12]);
	printf ("Forecast time in units : %d\n", readInt4(buf+13));
	break;
	}

	return 0;
	}

	int GRIB2_readSection5 (int fd, T_grib2Record *rec2, int secLen) {
	if (read (fd, buf, secLen) != secLen)
	return -1;

	printf ("==== Data Representation Section ====\n");
	printf ("Data Representation Template = %d (0:Grid Point Data - Simple Packing)\n", readInt2(buf+4));
	switch (readInt2(buf+4)) {
	case 0: // Grid Point Data - Simple Packing
	rec2->s5.R = ieee2flt(buf+6);
	rec2->s5.E = readSignedInt2(buf+10);
	rec2->s5.Epow2 = pow(2.,rec2->s5.E);
	rec2->s5.D = readSignedInt2(buf+12);
	rec2->s5.Dpow10 = pow(10.,rec2->s5.D);
	rec2->s5.nbBits = buf[14];
	printf ("Reference Value (R) = %f\n", rec2->s5.R);
	printf ("Binary Scale Factor (E) = %d\n", rec2->s5.E);
	printf ("Epow2 = %f\n", rec2->s5.Epow2);
	printf ("Decimal Scale Factor (D) = %d\n", rec2->s5.D);
	printf ("Dpow10 = %f\n", rec2->s5.Dpow10);
	printf ("nbBits : %d\n", rec2->s5.nbBits);
	break;
	}
	return 0;
	}

	int GRIB2_readSection7(int fd, T_grib2Record *rec2, int secLen) {
	if (read (fd, buf, secLen) != secLen)
	return -1;

	printf ("==== DATA SECTION ====\n");

	// Read data in the order given by isAdjacentI
	int i, j, x;
	double *pt = rec2->s7.data;
	unsigned int startbit = 0;

	for (j=0; j<rec2->s3.Nj; j++) {
	for (i=0; i<rec2->s3.Ni; i++) {
	x = readPackedBits(buf, startbit, rec2->s5.nbBits);
	pt = (rec2->s5.R + xrec2->s5.Epow2)/rec2->s5.Dpow10;
	printf ("data[%d,%d] = %d => %f\n", i, j, x, *pt);
	startbit += rec2->s5.nbBits;
	pt ++;
	}
	}
	return 0;
	}

	int GRIB2(int fd, T_grib2Record *rec2) {
	off_t startOfSection = lseek(fd, 0, SEEK_CUR);

	while (1) {
	if (read (fd, buf, 5) != 5)
	break;
	if (!memcmp(buf, "7777", 4)) {
	printf ("\nSection 8, 7777 End of Record\n");
	break;
	}
	long secLen = readInt4(buf);
	int secNum = buf[4];
	printf ("\nSection %d, len=%ld\n", secNum, secLen);

	switch (secNum) {
	case 1 :
	GRIB2_readSection1 (fd, rec2, secLen-5);
	break;
	case 3 :
	GRIB2_readSection3 (fd, rec2, secLen-5);
	break;
	case 4 :
	GRIB2_readSection4 (fd, rec2, secLen-5);
	break;
	case 5 :
	GRIB2_readSection5 (fd, rec2, secLen-5);
	break;
	case 7 :
	GRIB2_readSection7 (fd, rec2, secLen-5);
	break;
	}

	startOfSection = lseek (fd, startOfSection + secLen, SEEK_SET);
	}
	return 0;
	}

	// IS : Indicator Section
	int GRIBx_readSection0(int fd, T_gribRecord rec, T_grib2Record rec2) {
	rec->s0.offsetInFile = lseek (fd, 0L, SEEK_CUR);
	if (read (fd, buf, 8) != 8)
	return -1;
	memcpy (rec->s0.magic, buf, 4);
	rec->s0.magic[4] = 0;
	if (strcmp(rec->s0.magic, "GRIB"))
	return -2;
	rec->s0.editionNumber = buf[7];

	switch (rec->s0.editionNumber) {
	case 1:
	rec->s0.sectionSize = 8;
	rec->s0.totalLength = readInt3(buf+4);
	break;
	case 2:
	rec2->s0.offsetInFile = lseek (fd, 0L, SEEK_CUR);
	rec->s0.sectionSize = 16;
	rec2->s0.sectionSize = 16;
	rec2->s0.editionNumber = buf[7];
	strcpy (rec2->s0.magic, rec->s0.magic);
	rec2->s0.discipline = buf[6];
	if (read (fd, buf, 8) != 8)
	return -3;
	rec2->s0.totalLength = readLong8(buf);
	rec->s0.totalLength = rec2->s0.totalLength;
	break;
	}
	return 0;
	}

	void GRIBx_dumpSection0(T_gribRecord rec, T_grib2Record rec2) {
	if (!Verbose) return;
	printf ("\nIS :\n");
	printf ("s0.magic = %s\n", rec->s0.magic);
	printf ("s0.editionNumber = %d\n", rec->s0.editionNumber);
	switch (rec->s0.editionNumber) {
	case 1:
	printf ("s0.offsetInFile = %ld\n", rec->s0.offsetInFile);
	printf ("s0.sectionSize = %ld\n", rec->s0.sectionSize);
	printf ("s0.totalLength = %ld\n", rec->s0.totalLength);
	break;
	case 2:
	printf ("s0.offsetInFile = %ld\n", rec2->s0.offsetInFile);
	printf ("s0.sectionSize = %ld\n", rec2->s0.sectionSize);
	printf ("s0.discipline = %d\n", rec2->s0.discipline);
	printf ("s0.totalLength = %ld\n", rec2->s0.totalLength);
	break;
	}
	printf ("\n");
	}

	int GetGRIB(char fn, uint8_t dataType, uint8_t levelType, struct tm tmRequested, double La, double Lo, double *res) {
	T_gribRecord rec;
	T_grib2Record rec2;

	if (Verbose)
	printf ("REQUEST %04d/%02d/%02d at %02dh UTC : datatype:%d level=%d. for [%f,%f]\n",
	tmRequested->tm_year+1900, tmRequested->tm_mon+1, tmRequested->tm_mday, tmRequested->tm_hour,
	dataType, levelType, La, Lo);

	int fd = open(fn, 0);
	if (fd < 0)
	return -1;

	bool end = false;
	uint32_t recnum = 0;

	while (!end) {
	if (GRIBx_readSection0 (fd, &rec, &rec2) < 0) {
	end = true;
	break;
	}
	if (rec.s0.editionNumber == 1) {
	time_t t0;
	struct tm tm;

	if (GRIB1_readSection1 (fd, &rec) < 0) {
	end = true;
	break;
	}

	memset (&tm, 0, sizeof(tm));
	tm.tm_year = rec.s1.year - 1900;
	tm.tm_mon = rec.s1.month - 1;
	tm.tm_mday = rec.s1.day;
	tm.tm_hour = rec.s1.hour;
	t0 = timegm(&tm);

	if (rec.s1.forcastTimeUnit == 0) { // min
	t0 += rec.s1.periodP1*60;
	}
	else if (rec.s1.forcastTimeUnit == 1) { // hour
	t0 += rec.s1.periodP1*3600;
	}
	else if (rec.s1.forcastTimeUnit == 2) { // day
	t0 += rec.s1.periodP1360024;
	}

	//*tm = gmtime (&t0);
	memcpy (&tm, gmtime(&t0), sizeof(tm));

	GRIB1_readSection2 (fd, &rec);

	if (Verbose)
	printf ("Record #%d contains %04d/%02d/%02d at %02dh UTC : datatype:%d level=%d from [%f,%f] to [%f,%f]\n",
	recnum,
	tm.tm_year+1900, tm.tm_mon+1, tm.tm_mday, tm.tm_hour, rec.s1.dataType, rec.s1.levelType,
	rec.s2.La1, rec.s2.Lo1, rec.s2.La2, rec.s2.Lo2);

	if (rec.s1.dataType == dataType
	&& rec.s1.levelType == levelType
	&& tm.tm_year == tmRequested->tm_year && tm.tm_mon == tmRequested->tm_mon && tm.tm_mday == tmRequested->tm_mday
	&& tm.tm_hour == tmRequested->tm_hour) {

	//GRIB1_readSection2 (fd, &rec);

	if (inRectangle(&rec, La, Lo)) {
	if (Verbose)
	printf ("Found\n");
	GRIB1_dumpSection1(&rec);
	GRIB1_dumpSection2(&rec);

	GRIB1_readSection3 (fd, &rec);
	GRIB1_dumpSection3(&rec);

	GRIB1_readSection4 (fd, &rec);
	//GRIB1_dumpSection4(&rec);

	GRIB1_readSection5 (fd, &rec);

	close (fd);
	*res = getDataAtValues (&rec, La, Lo);
	return 0;
	}
	}

	if (GRIB1_seekToNextRecord(fd, &rec) < 0) {
	end = true;
	break;
	}
	}
	recnum ++;
	}
	close (fd);
	if (Verbose)
	printf ("-------------------- Not found -------------------\n");
	return -1;
	}

	void DumpGrib(char *fn) {
	T_gribRecord rec;
	T_grib2Record rec2;
	int recnum;

	int fd = open(fn, 0);
	if (fd < 0)
	exit (1);

	printf ("File %s\n", fn);
	recnum = 1;
	while (1) {
	if (GRIBx_readSection0 (fd, &rec, &rec2) < 0) {
	printf ("End of file\n");
	break;
	}
	printf ("Record #%d\n", recnum);
	GRIBx_dumpSection0 (&rec, &rec2);

	switch (rec.s0.editionNumber) {
	case 1:
	GRIB1 (fd, &rec);
	break;
	case 2:
	GRIB2 (fd, &rec2);
	break;
	}
	recnum ++;
	if (GRIB1_seekToNextRecord(fd, &rec) < 0)
	break;
	}

	close (fd);
	}

	// 2016-08-26T08:55:47
	int iso8601_to_tm(char iso, struct tm tm) {

	if (!iso \|\| strlen(iso) < 19)
	return -1;
	if ((iso[4] != '-') \|\| (iso[7] != '-') \|\| (iso[10] != 'T') \|\| (iso[13] != ':') \|\| (iso[16] != ':'))
	return -1;

	memset (tm, 0, sizeof(struct tm));

	tm->tm_year = (iso[0]-'0')1000 + (iso[1]-'0')100 + (iso[2]-'0') * 10 + (iso[3]-'0') - 1900;
	tm->tm_mon = (iso[5]-'0')*10 + (iso[6]-'0') - 1;
	tm->tm_mday = (iso[8]-'0')*10 + (iso[9]-'0');
	tm->tm_isdst = -1;
	tm->tm_hour = (iso[11]-'0')*10 + (iso[12]-'0');
	tm->tm_min = (iso[14]-'0')*10 + (iso[15]-'0');

	// Seconds are decimal
	char tmp[7];
	memcpy (tmp, iso+17, 6);
	tmp[6] = 0;
	double seconds = atof(tmp);
	tm->tm_sec = (int)seconds;

	return 0;
	}


	void demo(char *file) {
	struct tm tmRequested;
	double La, Lo;
	double uwind, vwind, temp;
	int rc;

	tmRequested.tm_year = 2016-1900;
	tmRequested.tm_mon = 9-1;
	tmRequested.tm_mday = 14;
	tmRequested.tm_hour = 6;
	La = 43;
	Lo = 5;
	Verbose = true;

	//DumpGrib (file);

	/*
	double pressure;
	rc = GetGRIB (file, DATATYPE_PRESSURE, LEVELTYPE_MSL, &tmRequested, La, Lo, &pressure);
	if (rc >= 0)
	printf ("===> Pressure : %f\n", pressure);
	*/

	rc = GetGRIB (file, DATATYPE_UWIND, LEVELTYPE_HEIGHT, &tmRequested, La, Lo, &uwind);
	if (rc >= 0) {
	printf ("===> uwind : %f m/s = %f Kts\n", uwind, uwind*1.94);

	rc = GetGRIB (file, DATATYPE_VWIND, LEVELTYPE_HEIGHT, &tmRequested, La, Lo, &vwind);
	printf ("===> vwind : %f m/s = %f Kts\n", vwind, vwind*1.94);

	double wind = sqrt(uwinduwind + vwindvwind);
	printf ("===> wind : %f m/s = %f Kts\n", wind, wind*1.94);

	double A = atan2(uwind, vwind) * 180 / M_PI;
	A = A + 180;
	printf (" Angle %f\n", A);
	}

	rc = GetGRIB (file, DATATYPE_TEMP, LEVELTYPE_MSL, &tmRequested, La, Lo, &temp);
	if (rc >= 0)
	printf ("===> temp : %fK, %f°C\n", temp, temp-273.15);

	rc = GetGRIB (file, DATATYPE_TEMP, LEVELTYPE_HEIGHT, &tmRequested, La, Lo, &temp);
	if (rc >= 0)
	printf ("===> temp : %fK, %f°C\n", temp, temp-273.15);
	}

	int main(int ac, char **av) {

	struct tm tmRequested;
	double La, Lo;
	int rc;

	if (ac < 2) exit (1);
	if (ac < 7) {
	demo(av[1]);
	return 0;
	}

	iso8601_to_tm(av[2], &tmRequested);
	La = atof(av[3]);
	Lo = atof(av[4]);
	uint8_t datatype = atoi(av[5]);
	uint8_t level = atoi(av[6]);

	Verbose = false;

	double value;
	rc = GetGRIB (av[1], datatype, level, &tmRequested, La, Lo, &value);
	if (rc >= 0)
	printf ("%f\n", value);
	else
	printf ("???\n");
	return 0;
	}

grib.cNo OneTemporaryActions

File Metadata

grib.cView Options

Event Timeline

grib.c
No OneTemporary
Actions

grib.c
View Options