Jpeg2000 data encoding and bitmap support

README.md

@@ -21,6 +21,8 @@ import { GRIB } from 'vgrib2'
 const gribFile = fs.readFileSync('./path/to/file.grib2')
 
 const grib = GRIB.parse(gribFile)
+// Lookup grid data point from the first packet
+const point = GRIB.lookupDataPoint(grib[0], -39.75, 146.15)
 // No Table Lookup (useful for adding new tables)
 const gribNoLookup = GRIB.parseNoLookup(gribFile)
 ```

package.json

@@ -27,16 +27,22 @@
   ],
   "name": "vgrib2",
   "author": "Alessandro Vecchi",
+  "contributors": [
+    "Dmitry Shirokov @runk"
+  ],
   "module": "dist/vgrib2.esm.js",
   "devDependencies": {
-    "@types/buffer-from": "^1.1.0",
     "tsdx": "^0.14.1",
     "tslib": "^2.3.0",
     "typescript": "^3.9.10"
   },
   "dependencies": {
-    "base64-js": "^1.3.1",
     "buffer": "^6.0.3",
-    "ieee754": "^1.2.1"
+    "jpeg2000": "^1.1.0"
+  },
+  "homepage": "https://github.com/veech/vgrib2",
+  "repository": {
+    "type": "git",
+    "url": "https://github.com/veech/vgrib2.git"
   }
 }

src/data/index.ts

@@ -2,22 +2,25 @@ import { Buffer } from 'buffer/'
 
 import { DataRepresentationSectionValues } from '../sections/section-5'
 
-import { simpleUnpacking } from './simple-unpacking'
+import { simpleUnpacking } from './simple'
+import { jpeg2000Unpacking } from './jpeg2000'
+import { BitMapSectionValues } from '../sections'
 
 /**
  * Converts data Buffer according to data representation section
  * @param drs Data Representation Section
  * @param data Data to be converted
  * @returns Converted data
  */
-export const convertData = (drs: DataRepresentationSectionValues, data: Buffer) => {
+export const convertData = (drs: DataRepresentationSectionValues, bms:BitMapSectionValues, data: Buffer) => {
   const { dataRepresentationTemplate } = drs.contents
 
   switch (dataRepresentationTemplate) {
     case 0:
       return simpleUnpacking(drs, data)
-
+    case 40:
+      return jpeg2000Unpacking(drs, bms, data)
     default:
       throw new Error(`Template 7.${dataRepresentationTemplate} not defined`)
   }
-}
+}

src/data/jpeg2000.ts

@@ -0,0 +1,57 @@
+import { Buffer } from 'buffer/'
+import { DataRepresentationSectionValues } from '../sections/section-5'
+import { JpxImage } from 'jpeg2000'
+import { BitMapSectionValues } from '../sections'
+
+/**
+ * Data Template 7.0
+ *
+ * [Read more...](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp7-40.shtml)
+ */
+export const jpeg2000Unpacking = (drs: DataRepresentationSectionValues, bms: BitMapSectionValues, data: Buffer) => {
+  const jpx = new JpxImage()
+  jpx.parse(data)
+
+  if (jpx.componentsCount !== 1) throw new Error('JPEG Decoder: Only single component is supported')
+  if (jpx.tiles.length !== 1) throw new Error('JPEG Decoder: Only single tile is supported')
+  if (jpx.tiles[0].height !== 1) throw new Error('JPEG Decoder: Only single row (1xN) is supported')
+
+  const { dataRepresentation } = drs.contents
+
+  const D = dataRepresentation.decimalScaleFactor
+  const R = dataRepresentation.referenceValue
+  const E = dataRepresentation.binaryScaleFactor
+
+  const DD = Math.pow(10, D)
+  const EE = Math.pow(2, E)
+
+  const result: (number | null)[] = []
+
+  // A bit map applies to this product
+  // [See more](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table6-0.shtml)
+  if (bms.contents.bitMapIndicator === 0) {
+    if (!bms.contents.bitMap) {
+      throw new Error('Bit map is not defined')
+    }
+
+    let k = 0;
+    for (const byte of bms.contents.bitMap) {
+      // Apply bit map to the data.
+      // Length of data values is often smaller than the bit map itself. Bitmap is used to 
+      // indicate which data values are present, 1 bit meaning is present, 0 bit meaning is missing.
+      // [Read more](https://confluence.ecmwf.int/display/UDOC/What+is+the+GRIB+bitmap+-+ecCodes+GRIB+FAQ)
+      for (let i = 0; i < 8; i++) {
+        if (byte & (1 << i)) {
+          result.push((R + jpx.tiles[0].items[k++] * EE) / DD)
+        } else {
+          result.push(null)
+        }
+      }
+    }
+  } else {
+    // Do not use `.map` on Uint8Array, as it clamps the values to 0-255
+    for (const byte of jpx.tiles[0].items) result.push((R + byte * EE) / DD)
+  }
+
+  return result
+}

src/index.ts

@@ -41,6 +41,26 @@ const parseNoLookup = (data: ArrayBuffer): Array<GRIBPacketValues> => {
   return packets
 }
 
-export const GRIB = { parse, parseNoLookup }
+/**
+ * Lookup data point from GRIB packet based on lat/lon coordinates.
+ *
+ * @param packet GRIB packet
+ * @param lat Latitude
+ * @param lon Longitude
+ */
+const lookupDataPoint = (packet: GRIBPacket, lat: number, lon: number): number | null => {
+  const { la1, lo1, la2, lo2, dx, dy, ny } = packet.gridDefinition
+  if (lat < la1 || lat > la2 || lon < lo1 || lon > lo2) return null
+
+  // lon (W-E) is x
+  // lat (N-S) is y
+  const x = Math.abs(Math.round((lo1 - lon) / dx))
+  const y = Math.abs(Math.round((la1 - lat) / dy))
+
+  const idx = y * ny + x
+  return packet.data[idx]
+}
+
+export const GRIB = { parse, parseNoLookup, lookupDataPoint }
 
 export * from './types/grib'