switch to equal area projection for rasters. Add toRasterProjection a…

…nd getRasterBoxSpherical helper methods

packages/geoprocessing/package.json

@@ -143,6 +143,7 @@
     "aws-sdk": "2.1399.0",
     "babel-jest": "^26.0.1",
     "babel-loader": "^8.2.2",
+    "bbox-fns": "^0.19.0",
     "bytes": "^3.1.0",
     "canonicalize": "^2.0.0",
     "chalk": "^4.1.1",
@@ -186,6 +187,7 @@
     "pbf": "^3.2.1",
     "polygon-clipping": "0.15.3",
     "pretty-bytes": "^5.3.0",
+    "proj4": "^2.9.2",
     "promptly": "^3.2.0",
     "rbush": "^3.0.1",
     "react": "^16.14.0",
@@ -196,6 +198,7 @@
     "react-table": "^7.7.0",
     "react-test-renderer": "^16.14.0",
     "read-pkg-up": "^7.0.1",
+    "reproject-geojson": "^0.5.0",
     "request": "^2.88.2",
     "serve-static": "^1.15.0",
     "slonik": "33.3.1",

packages/geoprocessing/scripts/base/datasources/importRasterDatasource.ts

@@ -68,7 +68,7 @@ export async function genCog(config: ImportRasterDatasourceConfig) {
   const warpDst = getCogPath(config.dstPath, config.datasourceId, "_4326");
   const dst = getCogPath(config.dstPath, config.datasourceId);
   // reproject
-  await $`gdalwarp -t_srs "EPSG:4326" --config GDAL_PAM_ENABLED NO --config GDAL_CACHEMAX 500 -wm 500 -multi -wo NUM_THREADS=ALL_CPUS ${src} ${warpDst}`;
+  await $`gdalwarp -t_srs "EPSG:6933" --config GDAL_PAM_ENABLED NO --config GDAL_CACHEMAX 500 -wm 500 -multi -wo NUM_THREADS=ALL_CPUS ${src} ${warpDst}`;
   // cloud-optimize
   await $`gdal_translate -b ${config.band} -r nearest --config GDAL_PAM_ENABLED NO --config GDAL_CACHEMAX 500 -co COMPRESS=LZW -co NUM_THREADS=ALL_CPUS -of COG -stats ${warpDst} ${dst}`;
   await $`rm ${warpDst}`;

packages/geoprocessing/scripts/base/datasources/precalcRasterDatasource.test.ts

@@ -127,15 +127,15 @@ describe("precalcRasterDatasource", () => {
 
     const sumMetric = firstMatchingMetric(metrics, (m) => m.metricId === "sum");
     expect(sumMetric).toBeTruthy();
-    expect(sumMetric.value).toBe(49);
+    expect(sumMetric.value).toBe(53);
 
     fs.removeSync(dsFilePath);
     fs.removeSync(path.join(dstPath, `${geogDatasourceId}.fgb`));
     fs.removeSync(path.join(dstPath, `${geogDatasourceId}.json`));
     fs.removeSync(geogFilePath);
     fs.removeSync(path.join(dstPath, `${datasourceId}.tif`));
     fs.removeSync(precalcFilePath);
-  }, 10000);
+  }, 60000);
 
   test("precalcRasterDatasource - multiple geog scenarios with external subdivided datasource", async () => {
     const dsFilename = "datasources_precalc_raster_test_9.json";
@@ -162,7 +162,7 @@ describe("precalcRasterDatasource", () => {
       },
     ]);
 
-    // Import raster to test against with known counts (49 Samoa and 17 American Samoa)
+    // Import raster to test against with known counts (both Samoa and American Samoa)
     await importDatasource(
       projectClient,
       {
@@ -254,19 +254,19 @@ describe("precalcRasterDatasource", () => {
       metrics,
       (m) => m.geographyId === "geog-box-filter"
     );
-    expect(boxFilterMetric.value).toEqual(91);
+    expect(boxFilterMetric.value).toEqual(101);
 
     const singleFilterMetric = firstMatchingMetric(
       metrics,
       (m) => m.geographyId === "geog-single-filter"
     );
-    expect(singleFilterMetric.value).toEqual(70);
+    expect(singleFilterMetric.value).toEqual(76);
 
     const doubleFilterMetric = firstMatchingMetric(
       metrics,
       (m) => m.geographyId === "geog-double-filter"
     );
-    expect(doubleFilterMetric.value).toEqual(88);
+    expect(doubleFilterMetric.value).toEqual(98);
 
     fs.removeSync(dsFilePath);
     fs.removeSync(path.join(dstPath, `${rasterDatasourceId}.tif`));
@@ -301,7 +301,7 @@ describe("precalcRasterDatasource", () => {
       },
     ]);
 
-    // Import raster to test against with known counts (49 Samoa and 17 American Samoa)
+    // Import raster to test against with known counts (both Samoa and American Samoa)
     await importDatasource(
       projectClient,
       {
@@ -398,19 +398,19 @@ describe("precalcRasterDatasource", () => {
       metrics,
       (m) => m.geographyId === "geog-box-filter"
     );
-    expect(boxFilterMetric.value).toEqual(69);
+    expect(boxFilterMetric.value).toEqual(70);
 
     const singleFilterMetric = firstMatchingMetric(
       metrics,
       (m) => m.geographyId === "geog-single-filter"
     );
-    expect(singleFilterMetric.value).toEqual(49);
+    expect(singleFilterMetric.value).toEqual(53);
 
     const doubleFilterMetric = firstMatchingMetric(
       metrics,
       (m) => m.geographyId === "geog-double-filter"
     );
-    expect(doubleFilterMetric.value).toEqual(66);
+    expect(doubleFilterMetric.value).toEqual(69);
 
     fs.removeSync(dsFilePath);
     fs.removeSync(path.join(dstPath, `${rasterDatasourceId}.tif`));
@@ -453,7 +453,7 @@ describe("precalcRasterDatasource", () => {
       }
     );
 
-    // Import raster to test against with known counts (49 Samoa and 21 American Samoa)
+    // Import raster to test against with known counts (both Samoa and American Samoa)
     await importDatasource(
       projectClient,
       {
@@ -545,19 +545,19 @@ describe("precalcRasterDatasource", () => {
       metrics,
       (m) => m.geographyId === "geog-box-filter"
     );
-    expect(noFilterMetric.value).toEqual(66);
+    expect(noFilterMetric.value).toEqual(69);
 
     const singleFilterMetric = firstMatchingMetric(
       metrics,
       (m) => m.geographyId === "geog-single-filter"
     );
-    expect(singleFilterMetric.value).toEqual(49);
+    expect(singleFilterMetric.value).toEqual(53);
 
     const doubleFilterMetric = firstMatchingMetric(
       metrics,
       (m) => m.geographyId === "geog-double-filter"
     );
-    expect(doubleFilterMetric.value).toEqual(66);
+    expect(doubleFilterMetric.value).toEqual(69);
 
     fs.removeSync(dsFilePath);
     fs.removeSync(path.join(dstPath, `${rasterDatasourceId}.tif`));

packages/geoprocessing/scripts/base/datasources/precalcRasterDatasource.ts

@@ -14,6 +14,7 @@ import {
   BBox,
   ProjectClientBase,
   datasourceConfig,
+  getRasterBoxSpherical,
 } from "../../../src";
 import bbox from "@turf/bbox";
 
@@ -83,7 +84,7 @@ export async function genRasterMetrics(
   const dstPath = extraOptions.newDstPath || datasourceConfig.defaultDstPath;
 
   // console.log(
-  //   `DATASOURCE: ${datasource.datasourceId}}, GEOGRAPHY: ${geography.geographyId}}\n`
+  //   `DATASOURCE: ${datasource.datasourceId}, GEOGRAPHY: ${geography.geographyId}\n`
   // );
 
   // Reads in geography vector data as FeatureCollection
@@ -95,7 +96,7 @@ export async function genRasterMetrics(
   );
   // console.log("geographyFeatureColl", JSON.stringify(geographyFeatureColl));
 
-  const rasterBbox: BBox = [raster.xmin, raster.ymin, raster.xmax, raster.ymax];
+  const rasterBbox: BBox = getRasterBoxSpherical(raster);
 
   // If there's no overlap between geography and raster, return empty metric
   if (!bboxOverlap(bbox(geographyFeatureColl), rasterBbox)) {

packages/geoprocessing/src/toolbox/geoblaze/geoblaze.ts

@@ -1,6 +1,9 @@
 import { Polygon, Histogram } from "../../types";
-import { Feature, MultiPolygon, FeatureCollection } from "@turf/helpers";
+import { Feature, MultiPolygon, FeatureCollection, BBox } from "@turf/helpers";
 import geoblaze, { Georaster } from "geoblaze";
+import reprojectGeoJSONPlugable from "reproject-geojson/pluggable.js";
+import proj4 from "proj4";
+import { reproject } from "bbox-fns";
 
 /**
  * Returns sum of value overlap with geometry.  If no cells with a value are found within the geometry overlap, returns 0.
@@ -12,8 +15,9 @@ export const getSum = async (
     | FeatureCollection<Polygon | MultiPolygon>
 ) => {
   let sum = 0;
+  const finalFeat = toRasterProjection(raster, feat);
   try {
-    const result = await geoblaze.sum(raster, feat);
+    const result = await geoblaze.sum(raster, finalFeat);
     sum = result[0];
   } catch (err) {
     console.log(
@@ -49,3 +53,61 @@ export const getHistogram = async (
   }
   return histogram;
 };
+
+/**
+ * Reprojects a feature to the same projection as the raster.
+ * @param raster
+ * @param feat
+ * @returns feature in projection of raster
+ * @throws if raster projection is not 4326 (backwards-compatibility) or 6933 (new equal area)
+ */
+export const toRasterProjection = (
+  raster: Georaster,
+  feat?:
+    | Feature<Polygon | MultiPolygon>
+    | FeatureCollection<Polygon | MultiPolygon>
+) => {
+  if (raster.projection === 4326) {
+    return feat;
+  } else if (raster.projection === 6933) {
+    const { forward } = proj4("EPSG:4326", "EPSG:6933");
+    proj4.defs(
+      "EPSG:6933",
+      "+proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"
+    );
+    return reprojectGeoJSONPlugable(feat, {
+      reproject: forward,
+    });
+  } else {
+    throw new Error(
+      `Unexpected projection for raster: ${raster.projection}. Expected 6933 or 4326`
+    );
+  }
+};
+
+/**
+ * @param raster
+ * @returns the bounding box of the raster in WGS84 lat/lon coordinates
+ */
+export const getRasterBoxSpherical = (raster: Georaster) => {
+  if (raster.projection === 4326) {
+    const bbox: BBox = [raster.xmin, raster.ymin, raster.xmax, raster.ymax];
+    return bbox;
+  } else if (raster.projection === 6933) {
+    // Reproject back to spherical coordinates
+    proj4.defs(
+      "EPSG:6933",
+      "+proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"
+    );
+    const { inverse } = proj4("EPSG:4326", "EPSG:6933");
+    const rasterBbox: BBox = reproject(
+      [raster.xmin, raster.ymin, raster.xmax, raster.ymax],
+      inverse
+    );
+    return rasterBbox;
+  } else {
+    throw new Error(
+      `Unexpected projection for raster: ${raster.projection}. Expected 6933 or 4326`
+    );
+  }
+};