From fb776645301e897a0ad3eff0a746e0be26de9d09 Mon Sep 17 00:00:00 2001
From: amelia <asnyder@usgs.gov>
Date: Tue, 6 Feb 2024 16:04:33 -0600
Subject: [PATCH] add WGS84 conversion for iclus_hc
---
catalog/iclus_hc/collection.json | 12 +--
...iclus_hc_create_collection_from_zarr.ipynb | 75 +++++++++++++------
2 files changed, 60 insertions(+), 27 deletions(-)
diff --git a/catalog/iclus_hc/collection.json b/catalog/iclus_hc/collection.json
index 70ae3401..5f5fce90 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\":\"ProjectedCRS\",\"name\":\"unknown\",\"base_crs\":{\"name\":\"unknown\",\"datum\":{\"type\":\"GeodeticReferenceFrame\",\"name\":\"unknown\",\"ellipsoid\":{\"name\":\"GRS 1980\",\"semi_major_axis\":6378137,\"inverse_flattening\":298.257222101}},\"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\":-96,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8822}},{\"name\":\"Latitude of 1st standard parallel\",\"value\":2,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8823}},{\"name\":\"Latitude of 2nd standard parallel\",\"value\":9,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8824}},{\"name\":\"Easting at false origin\",\"value\":0,\"unit\":\"metre\",\"id\":{\"authority\":\"EPSG\",\"code\":8826}},{\"name\":\"Northing at false origin\",\"value\":0,\"unit\":\"metre\",\"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\"}]}}"
+ "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\":[]}"
},
"y": {
"type": "spatial",
@@ -39,7 +39,7 @@
3172445.2073140107
],
"step": -100.0,
- "reference_system": "{\"$schema\":\"https://proj.org/schemas/v0.5/projjson.schema.json\",\"type\":\"ProjectedCRS\",\"name\":\"unknown\",\"base_crs\":{\"name\":\"unknown\",\"datum\":{\"type\":\"GeodeticReferenceFrame\",\"name\":\"unknown\",\"ellipsoid\":{\"name\":\"GRS 1980\",\"semi_major_axis\":6378137,\"inverse_flattening\":298.257222101}},\"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\":-96,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8822}},{\"name\":\"Latitude of 1st standard parallel\",\"value\":2,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8823}},{\"name\":\"Latitude of 2nd standard parallel\",\"value\":9,\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8824}},{\"name\":\"Easting at false origin\",\"value\":0,\"unit\":\"metre\",\"id\":{\"authority\":\"EPSG\",\"code\":8826}},{\"name\":\"Northing at false origin\",\"value\":0,\"unit\":\"metre\",\"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\"}]}}"
+ "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\":[]}"
}
},
"cube:variables": {
@@ -909,10 +909,10 @@
"spatial": {
"bbox": [
[
- -2356064.09086355,
- 269345.2073140107,
- 2258235.90913645,
- 3172445.2073140107
+ -118.73555823136944,
+ 22.872641668268155,
+ -65.34617124487134,
+ 48.24390993761895
]
]
},
diff --git a/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb b/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb
index 8abe442f..4f023f4e 100644
--- a/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb
+++ b/workflows/archive/iclus_hc_create_collection_from_zarr.ipynb
@@ -38,6 +38,7 @@
"import cfunits\n",
"import json\n",
"import sys\n",
+ "from pyproj import Transformer\n",
"sys.path.insert(1, '..')\n",
"import stac_helpers"
]
@@ -214,7 +215,28 @@
"outputs": [],
"source": [
"ds = ds.metpy.parse_cf()\n",
- "crs = ds[list(ds.keys())[0]].metpy.cartopy_crs"
+ "crs = ds[list(ds.keys())[0]].metpy.pyproj_proj\n",
+ "cartopycrs = ds[list(ds.keys())[0]].metpy.cartopy_crs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "391b86aa-c5a8-4bae-9303-3a40ede8e8b6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds.crs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6bc00354-50d0-4e4c-ae42-1964fe7acba5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "crs.to_proj4()"
]
},
{
@@ -246,13 +268,36 @@
"# left, bottom, right, top\n",
"\n",
"# Note: try changing around the commented out lines below to get type float ratherthan a numpy float\n",
- "#coord_bounds = [ds[dim_names_dict['X']].data.min().compute().astype(float), ds[dim_names_dict['Y']].data.min().compute().astype(float), ds[dim_names_dict['X']].data.max().compute().astype(float), ds[dim_names_dict['Y']].data.max().compute().astype(float)]\n",
- "#coord_bounds = [ds[dim_names_dict['X']].data.min().compute().astype(float).tolist(), ds[dim_names_dict['Y']].data.min().compute().astype(float).tolist(), ds[dim_names_dict['X']].data.max().compute().astype(float).tolist(), ds[dim_names_dict['Y']].data.max().compute().astype(float).tolist()]\n",
- "coord_bounds = [ds[dim_names_dict['X']].data.min().astype(float).item(), ds[dim_names_dict['Y']].data.min().astype(float).item(), ds[dim_names_dict['X']].data.max().astype(float).item(), ds[dim_names_dict['Y']].data.max().astype(float).item()]\n",
- "print(coord_bounds)\n",
- "print(f'\\ncoord_bounds data type: {type(coord_bounds[0])}')\n",
+ "#spatial_bounds = [ds[dim_names_dict['X']].data.min().compute().astype(float), ds[dim_names_dict['Y']].data.min().compute().astype(float), ds[dim_names_dict['X']].data.max().compute().astype(float), ds[dim_names_dict['Y']].data.max().compute().astype(float)]\n",
+ "#spatial_bounds = [ds[dim_names_dict['X']].data.min().compute().astype(float).tolist(), ds[dim_names_dict['Y']].data.min().compute().astype(float).tolist(), ds[dim_names_dict['X']].data.max().compute().astype(float).tolist(), ds[dim_names_dict['Y']].data.max().compute().astype(float).tolist()]\n",
+ "spatial_bounds = [ds[dim_names_dict['X']].data.min().astype(float).item(), ds[dim_names_dict['Y']].data.min().astype(float).item(), ds[dim_names_dict['X']].data.max().astype(float).item(), ds[dim_names_dict['Y']].data.max().astype(float).item()]\n",
+ "print(spatial_bounds)\n",
+ "print(f'\\ncoord_bounds data type: {type(spatial_bounds[0])}')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1426d661-40b8-469a-a3d8-74f926e1c30c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "transformer = Transformer.from_crs(crs.crs, \"EPSG:4326\", always_xy=True)\n",
+ "lon1, lat1 = transformer.transform(xx=spatial_bounds[0], yy=spatial_bounds[1])\n",
+ "lon2, lat2 = transformer.transform(xx=spatial_bounds[2], yy=spatial_bounds[3])\n",
+ "print(f'lower left coordinates (WGS84): {lon1}, {lat1}')\n",
+ "print(f'upper right coordinates (WGS84): {lon2}, {lat2}')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "72eb51cc-aff2-4b53-8ee2-df183a59d3d5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
"# create a spatial extent object \n",
- "spatial_extent = pystac.SpatialExtent(bboxes=[coord_bounds])"
+ "spatial_extent = pystac.SpatialExtent(bboxes=[[lon1, lat1, lon2, lat2]])"
]
},
{
@@ -418,18 +463,6 @@
"print(dims)"
]
},
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "00a18a29-fb9a-4b56-8009-493122997b16",
- "metadata": {},
- "outputs": [],
- "source": [
- "# get x, y bounds for extent of those dimensions (required)\n",
- "xy_bounds = [ds[dim_names_dict['X']].data.min().astype(float).item(), ds[dim_names_dict['Y']].data.min().astype(float).item(), ds[dim_names_dict['X']].data.max().astype(float).item(), ds[dim_names_dict['Y']].data.max().astype(float).item()]\n",
- "print(xy_bounds)"
- ]
- },
{
"cell_type": "markdown",
"id": "e7dc357c-91ec-49ae-83e5-400f791f9792",
@@ -656,8 +689,8 @@
"\n",
"# we do not recommend including redundant dimensions (do not include x,y if you have lon,lat)\n",
"# note that the extent of each dimension should be pulled from the dataset\n",
- "dims_dict = {dim_names_dict['X']: pystac.extensions.datacube.Dimension({'type': 'spatial', 'axis': 'x', 'description': stac_helpers.get_long_name(ds, dim_names_dict['X']), 'extent': [xy_bounds[0], xy_bounds[2]], 'step': x_step, 'reference_system': projjson}),\n",
- " dim_names_dict['Y']: pystac.extensions.datacube.Dimension({'type': 'spatial', 'axis': 'y', 'description': stac_helpers.get_long_name(ds, dim_names_dict['Y']), 'extent': [xy_bounds[1], xy_bounds[3]], 'step': y_step, 'reference_system': projjson}),\n",
+ "dims_dict = {dim_names_dict['X']: pystac.extensions.datacube.Dimension({'type': 'spatial', 'axis': 'x', 'description': stac_helpers.get_long_name(ds, dim_names_dict['X']), 'extent': [spatial_bounds[0], spatial_bounds[2]], 'step': x_step, 'reference_system': projjson}),\n",
+ " dim_names_dict['Y']: pystac.extensions.datacube.Dimension({'type': 'spatial', 'axis': 'y', 'description': stac_helpers.get_long_name(ds, dim_names_dict['Y']), 'extent': [spatial_bounds[1], spatial_bounds[3]], 'step': y_step, 'reference_system': projjson}),\n",
"}"
]
},
--
GitLab