1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
|
#include <stdlib.h>
#include <math.h>
#include "grib2.h"
void simpack(g2float *fld,g2int ndpts,g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
//$$$ SUBPROGRAM DOCUMENTATION BLOCK
// . . . .
// SUBPROGRAM: simpack
// PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-11-06
//
// ABSTRACT: This subroutine packs up a data field using the simple
// packing algorithm as defined in the GRIB2 documention. It
// also fills in GRIB2 Data Representation Template 5.0 with the
// appropriate values.
//
// PROGRAM HISTORY LOG:
// 2002-11-06 Gilbert
//
// USAGE: CALL simpack(fld,ndpts,idrstmpl,cpack,lcpack)
// INPUT ARGUMENT LIST:
// fld[] - Contains the data values to pack
// ndpts - The number of data values in array fld[]
// idrstmpl - Contains the array of values for Data Representation
// Template 5.0
// [0] = Reference value - ignored on input
// [1] = Binary Scale Factor
// [2] = Decimal Scale Factor
// [3] = Number of bits used to pack data, if value is
// > 0 and <= 31.
// If this input value is 0 or outside above range
// then the num of bits is calculated based on given
// data and scale factors.
// [4] = Original field type - currently ignored on input
// Data values assumed to be reals.
//
// OUTPUT ARGUMENT LIST:
// idrstmpl - Contains the array of values for Data Representation
// Template 5.0
// [0] = Reference value - set by simpack routine.
// [1] = Binary Scale Factor - unchanged from input
// [2] = Decimal Scale Factor - unchanged from input
// [3] = Number of bits used to pack data, unchanged from
// input if value is between 0 and 31.
// If this input value is 0 or outside above range
// then the num of bits is calculated based on given
// data and scale factors.
// [4] = Original field type - currently set = 0 on output.
// Data values assumed to be reals.
// cpack - The packed data field
// lcpack - length of packed field starting at cpack.
//
// REMARKS: None
//
// ATTRIBUTES:
// LANGUAGE: C
// MACHINE:
//
//$$$
{
static g2int zero=0;
g2int *ifld;
g2int j,nbits,imin,imax,maxdif,nbittot,left;
g2float bscale,dscale,rmax,rmin,temp;
double maxnum;
static g2float alog2=0.69314718; // ln(2.0)
bscale=int_power(2.0,-idrstmpl[1]);
dscale=int_power(10.0,idrstmpl[2]);
if (idrstmpl[3] <= 0 || idrstmpl[3] > 31)
nbits=0;
else
nbits=idrstmpl[3];
//
// Find max and min values in the data
//
rmax=fld[0];
rmin=fld[0];
for (j=1;j<ndpts;j++) {
if (fld[j] > rmax) rmax=fld[j];
if (fld[j] < rmin) rmin=fld[j];
}
ifld=calloc(ndpts,sizeof(g2int));
//
// If max and min values are not equal, pack up field.
// If they are equal, we have a constant field, and the reference
// value (rmin) is the value for each point in the field and
// set nbits to 0.
//
if (rmin != rmax) {
//
// Determine which algorithm to use based on user-supplied
// binary scale factor and number of bits.
//
if (nbits==0 && idrstmpl[1]==0) {
//
// No binary scaling and calculate minumum number of
// bits in which the data will fit.
//
imin=(g2int)rint(rmin*dscale);
imax=(g2int)rint(rmax*dscale);
maxdif=imax-imin;
temp=log((double)(maxdif+1))/alog2;
nbits=(g2int)ceil(temp);
rmin=(g2float)imin;
// scale data
for(j=0;j<ndpts;j++)
ifld[j]=(g2int)rint(fld[j]*dscale)-imin;
}
else if (nbits!=0 && idrstmpl[1]==0) {
//
// Use minimum number of bits specified by user and
// adjust binary scaling factor to accomodate data.
//
rmin=rmin*dscale;
rmax=rmax*dscale;
maxnum=int_power(2.0,nbits)-1;
temp=log(maxnum/(rmax-rmin))/alog2;
idrstmpl[1]=(g2int)ceil(-1.0*temp);
bscale=int_power(2.0,-idrstmpl[1]);
// scale data
for (j=0;j<ndpts;j++)
ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
}
else if (nbits==0 && idrstmpl[1]!=0) {
//
// Use binary scaling factor and calculate minumum number of
// bits in which the data will fit.
//
rmin=rmin*dscale;
rmax=rmax*dscale;
maxdif=(g2int)rint((rmax-rmin)*bscale);
temp=log((double)(maxdif+1))/alog2;
nbits=(g2int)ceil(temp);
// scale data
for (j=0;j<ndpts;j++)
ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
}
else if (nbits!=0 && idrstmpl[1]!=0) {
//
// Use binary scaling factor and use minumum number of
// bits specified by user. Dangerous - may loose
// information if binary scale factor and nbits not set
// properly by user.
//
rmin=rmin*dscale;
// scale data
for (j=0;j<ndpts;j++)
ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
}
//
// Pack data, Pad last octet with Zeros, if necessary,
// and calculate the length of the packed data in bytes
//
sbits(cpack,ifld+0,0,nbits,0,ndpts);
nbittot=nbits*ndpts;
left=8-(nbittot%8);
if (left != 8) {
sbit(cpack,&zero,nbittot,left); // Pad with zeros to fill Octet
nbittot=nbittot+left;
}
*lcpack=nbittot/8;
}
else {
nbits=0;
*lcpack=0;
}
//
// Fill in ref value and number of bits in Template 5.0
//
//printf("SAGmkieee %f\n",rmin);
mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
//printf("SAGmkieee %ld\n",idrstmpl[0]);
idrstmpl[3]=nbits;
idrstmpl[4]=0; // original data were reals
free(ifld);
}
|