From 1da985c846221626bef80e9ebbde8379f28961b0 Mon Sep 17 00:00:00 2001
From: amsnyder <asnyder@usgs.gov>
Date: Mon, 8 Apr 2024 13:00:34 -0500
Subject: [PATCH] fix crs var detection

---
 catalog/iclus_hc/collection.json              |  6 +--
 ...iclus_hc_create_collection_from_zarr.ipynb | 41 ++++++++++++-------
 2 files changed, 30 insertions(+), 17 deletions(-)

diff --git a/catalog/iclus_hc/collection.json b/catalog/iclus_hc/collection.json
index 721bd8e1..66a79f91 100644
--- a/catalog/iclus_hc/collection.json
+++ b/catalog/iclus_hc/collection.json
@@ -28,7 +28,7 @@
         2258235.90913645
       ],
       "step": 100.0,
-      "reference_system": "{\"$schema\":\"https://proj.org/schemas/v0.5/projjson.schema.json\",\"type\":\"Conversion\",\"name\":\"PROJ-based coordinate operation\",\"method\":{\"name\":\"PROJ-based operation method: +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs\"},\"parameters\":[]}"
+      "reference_system": "{\"$schema\":\"https://proj.org/schemas/v0.5/projjson.schema.json\",\"type\":\"ProjectedCRS\",\"name\":\"undefined\",\"base_crs\":{\"name\":\"undefined\",\"datum\":{\"type\":\"GeodeticReferenceFrame\",\"name\":\"undefined\",\"ellipsoid\":{\"name\":\"undefined\",\"semi_major_axis\":6378137,\"inverse_flattening\":298.257222101},\"prime_meridian\":{\"name\":\"undefined\",\"longitude\":0}},\"coordinate_system\":{\"subtype\":\"ellipsoidal\",\"axis\":[{\"name\":\"Longitude\",\"abbreviation\":\"lon\",\"direction\":\"east\",\"unit\":\"degree\"},{\"name\":\"Latitude\",\"abbreviation\":\"lat\",\"direction\":\"north\",\"unit\":\"degree\"}]}},\"conversion\":{\"name\":\"unknown\",\"method\":{\"name\":\"Albers Equal Area\",\"id\":{\"authority\":\"EPSG\",\"code\":9822}},\"parameters\":[{\"name\":\"Latitude of false origin\",\"value\":23,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8821}},{\"name\":\"Longitude of false origin\",\"value\":-95.9999999999999,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8822}},{\"name\":\"Latitude of 1st standard parallel\",\"value\":29.5,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8823}},{\"name\":\"Latitude of 2nd standard parallel\",\"value\":45.5,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8824}},{\"name\":\"Easting at false origin\",\"value\":0,\"unit\":{\"type\":\"LinearUnit\",\"name\":\"Metre\",\"conversion_factor\":1},\"id\":{\"authority\":\"EPSG\",\"code\":8826}},{\"name\":\"Northing at false origin\",\"value\":0,\"unit\":{\"type\":\"LinearUnit\",\"name\":\"Metre\",\"conversion_factor\":1},\"id\":{\"authority\":\"EPSG\",\"code\":8827}}]},\"coordinate_system\":{\"subtype\":\"Cartesian\",\"axis\":[{\"name\":\"Easting\",\"abbreviation\":\"E\",\"direction\":\"east\",\"unit\":\"metre\"},{\"name\":\"Northing\",\"abbreviation\":\"N\",\"direction\":\"north\",\"unit\":\"metre\"}]}}"
     },
     "y": {
       "type": "spatial",
@@ -39,13 +39,13 @@
         3172445.2073140107
       ],
       "step": -100.0,
-      "reference_system": "{\"$schema\":\"https://proj.org/schemas/v0.5/projjson.schema.json\",\"type\":\"Conversion\",\"name\":\"PROJ-based coordinate operation\",\"method\":{\"name\":\"PROJ-based operation method: +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs\"},\"parameters\":[]}"
+      "reference_system": "{\"$schema\":\"https://proj.org/schemas/v0.5/projjson.schema.json\",\"type\":\"ProjectedCRS\",\"name\":\"undefined\",\"base_crs\":{\"name\":\"undefined\",\"datum\":{\"type\":\"GeodeticReferenceFrame\",\"name\":\"undefined\",\"ellipsoid\":{\"name\":\"undefined\",\"semi_major_axis\":6378137,\"inverse_flattening\":298.257222101},\"prime_meridian\":{\"name\":\"undefined\",\"longitude\":0}},\"coordinate_system\":{\"subtype\":\"ellipsoidal\",\"axis\":[{\"name\":\"Longitude\",\"abbreviation\":\"lon\",\"direction\":\"east\",\"unit\":\"degree\"},{\"name\":\"Latitude\",\"abbreviation\":\"lat\",\"direction\":\"north\",\"unit\":\"degree\"}]}},\"conversion\":{\"name\":\"unknown\",\"method\":{\"name\":\"Albers Equal Area\",\"id\":{\"authority\":\"EPSG\",\"code\":9822}},\"parameters\":[{\"name\":\"Latitude of false origin\",\"value\":23,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8821}},{\"name\":\"Longitude of false origin\",\"value\":-95.9999999999999,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8822}},{\"name\":\"Latitude of 1st standard parallel\",\"value\":29.5,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8823}},{\"name\":\"Latitude of 2nd standard parallel\",\"value\":45.5,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8824}},{\"name\":\"Easting at false origin\",\"value\":0,\"unit\":{\"type\":\"LinearUnit\",\"name\":\"Metre\",\"conversion_factor\":1},\"id\":{\"authority\":\"EPSG\",\"code\":8826}},{\"name\":\"Northing at false origin\",\"value\":0,\"unit\":{\"type\":\"LinearUnit\",\"name\":\"Metre\",\"conversion_factor\":1},\"id\":{\"authority\":\"EPSG\",\"code\":8827}}]},\"coordinate_system\":{\"subtype\":\"Cartesian\",\"axis\":[{\"name\":\"Easting\",\"abbreviation\":\"E\",\"direction\":\"east\",\"unit\":\"metre\"},{\"name\":\"Northing\",\"abbreviation\":\"N\",\"direction\":\"north\",\"unit\":\"metre\"}]}}"
     }
   },
   "cube:variables": {
     "crs": {
       "dimensions": [],
-      "type": "data",
+      "type": "auxiliary",
       "description": null,
       "unit": null
     },
diff --git a/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb b/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb
index 47b207f2..1b5e9e32 100644
--- a/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb
+++ b/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb
@@ -35,17 +35,14 @@
     "import pandas as pd\n",
     "import json\n",
     "import numpy as np\n",
-    "import metpy\n",
+    "import pyproj\n",
+    "from pyproj import Transformer\n",
     "import cartopy.crs as ccrs\n",
     "import cfunits\n",
     "import json\n",
     "import sys\n",
-    "from pyproj import Transformer\n",
     "sys.path.insert(1, '..')\n",
-    "import stac_helpers\n",
-    "import dask\n",
-    "import dask.array as da\n",
-    "from dask.distributed import Client, LocalCluster"
+    "import stac_helpers"
    ]
   },
   {
@@ -218,6 +215,18 @@
     "## Get crs info"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b2520710-a8a0-466d-adef-9ed01ed3dccf",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "crs_var = 'crs'"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -227,9 +236,13 @@
    },
    "outputs": [],
    "source": [
-    "ds  = ds.metpy.parse_cf()\n",
-    "crs = ds[list(ds.keys())[0]].metpy.pyproj_proj\n",
-    "#cartopycrs = ds[list(ds.keys())[0]].metpy.cartopy_crs"
+    "# use pyproj to automatically extract crs info\n",
+    "crs = pyproj.CRS.from_cf(ds[crs_var].attrs)\n",
+    "\n",
+    "# alternatively, create the appropriate cartopy projection\n",
+    "# crs = ccrs.LambertConformal(central_longitude=crs_info.longitude_of_central_meridian, \n",
+    "#                             central_latitude=crs_info.latitude_of_projection_origin,\n",
+    "#                             standard_parallels=crs_info.standard_parallel)"
    ]
   },
   {
@@ -249,7 +262,7 @@
    },
    "outputs": [],
    "source": [
-    "ds.crs"
+    "ds[crs_var]"
    ]
   },
   {
@@ -311,8 +324,8 @@
    },
    "outputs": [],
    "source": [
-    "# this is currently crashing - working on a dask application of this\n",
-    "XX, YY = dask.array.meshgrid(ds[dim_names_dict['X']].data, ds[dim_names_dict['Y']].data)"
+    "# uncomment if you wish to use dask\n",
+    "# XX, YY = dask.array.meshgrid(ds[dim_names_dict['X']].data, ds[dim_names_dict['Y']].data)"
    ]
   },
   {
@@ -342,7 +355,7 @@
    },
    "outputs": [],
    "source": [
-    "transformer = Transformer.from_crs(crs.crs, \"EPSG:4326\", always_xy=True)"
+    "transformer = Transformer.from_crs(crs, \"EPSG:4326\", always_xy=True)"
    ]
   },
   {
@@ -898,7 +911,7 @@
     "vars_dict={}\n",
     "for v in vars:\n",
     "    unit = stac_helpers.get_unit(ds, v)\n",
-    "    var_type = stac_helpers.get_var_type(ds, v)\n",
+    "    var_type = stac_helpers.get_var_type(ds, v, crs_var)\n",
     "    long_name = stac_helpers.get_long_name(ds, v)\n",
     "    vars_dict[v] = pystac.extensions.datacube.Variable({'dimensions':list(ds[v].dims), 'type': var_type, 'description': long_name, 'unit': unit})"
    ]
-- 
GitLab