From 2b415508a5eb6fdaf0f9d090abc977e662cfca7c Mon Sep 17 00:00:00 2001 From: spestana Date: Thu, 26 Sep 2024 20:46:45 +0000 Subject: [PATCH 1/3] Further refactoring --- src/goes_ortho/clip.py | 13 ++++++++++++- src/goes_ortho/get_data.py | 39 +++++++++++++++++++++++++++++++++++++- 2 files changed, 50 insertions(+), 2 deletions(-) diff --git a/src/goes_ortho/clip.py b/src/goes_ortho/clip.py index b6b0e6c..1ac7ad2 100644 --- a/src/goes_ortho/clip.py +++ b/src/goes_ortho/clip.py @@ -2,6 +2,7 @@ Functions for clipping GOES ABI imagery to smaller areas """ +import datetime as dt import logging import xarray as xr @@ -11,7 +12,7 @@ def subsetNetCDF(filepath, bounds, newfilepath=None): """ - Crop a GOES-R ABI NetCDF file to latitude/longitude bounds. + Crop a GOES-R ABI NetCDF file to latitude/longitude bounds, add datetime dims/coords. Parameters ------------ @@ -85,6 +86,16 @@ def subsetNetCDF(filepath, bounds, newfilepath=None): # Use these coordinates to subset the whole dataset y_rad_bnds, x_rad_bnds = [y_rad_n, y_rad_s], [x_rad_w, x_rad_e] ds = f.sel(x=slice(*x_rad_bnds), y=slice(*y_rad_bnds)) + + # while we have the file open, add datetime dims/coords to the file + # parse the start time from the file name (the part "s2022__________") + this_datetime = dt.datetime.strptime( + filepath.stem.split("_")[3][1:-1], "%Y%j%H%M%S" + ) + ds = ds.assign_coords({"time": this_datetime}) + ds = ds.expand_dims("time") + ds = ds.reset_coords(drop=True) + # Close the original file f.close() if newfilepath is None: diff --git a/src/goes_ortho/get_data.py b/src/goes_ortho/get_data.py index b28a330..14a48d5 100644 --- a/src/goes_ortho/get_data.py +++ b/src/goes_ortho/get_data.py @@ -130,7 +130,7 @@ def build_zarr(download_request_filepath, downloader="goes2go"): del source_group del source_array print("Done.") - return None + return outputFilepath def make_request_json( @@ -651,3 +651,40 @@ def bounds_from_geojson(geojson_filepath: str) -> Union[List[List[float]], List[ ] # if there is only one set of bounds, remove it from the list return bounds + + +def multi_nc_to_zarr( + nc_filepaths: List[Union[Path, str]], zarr_filepath: str +) -> Union[Path, str]: + """ + Read a list of GOES-R ABI netcdf filepaths, merge into a single zarr file along time dimension. + """ + + with go.io.dask_start_cluster( + workers=6, + threads=2, + open_browser=False, + verbose=True, + ) as _: + # read all the netcdf files, drop attributes where there are conflicts + ds = xr.open_mfdataset( + nc_filepaths, chunks={"time": 500}, combine_attrs="drop_conflicts" + ) + + # configure chunking on variables with dimensions of (time, x, y) + for variable, _ in ds.data_vars.items(): + if ds[variable].dims == ("time", "y", "x"): + ds[variable].data.rechunk( + {0: -1, 1: "auto", 2: "auto"}, block_size_limit=1e8, balance=True + ) + # Assign the dimensions of a chunk to variables to use for encoding afterwards + t, y, x = ( + ds[variable].data.chunks[0][0], + ds[variable].data.chunks[1][0], + ds[variable].data.chunks[2][0], + ) + ds[variable].encoding = {"chunks": (t, y, x)} + + ds.to_zarr(zarr_filepath) + + return zarr_filepath From 2b7358030a8b6b6744cceefb94c76c6e1221198f Mon Sep 17 00:00:00 2001 From: spestana Date: Fri, 27 Sep 2024 02:31:47 +0000 Subject: [PATCH 2/3] Update functions for working with zarr files, demo notebook --- docs/examples/example-workflow.ipynb | 534 +++++++++++++++++++++++++++ src/goes_ortho/get_data.py | 11 +- src/goes_ortho/orthorectify.py | 347 +++++++++++++++++ 3 files changed, 891 insertions(+), 1 deletion(-) create mode 100644 docs/examples/example-workflow.ipynb diff --git a/docs/examples/example-workflow.ipynb b/docs/examples/example-workflow.ipynb new file mode 100644 index 0000000..4001c31 --- /dev/null +++ b/docs/examples/example-workflow.ipynb @@ -0,0 +1,534 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "3aea64b8-4c4c-478f-85a3-a5c36321c4c9", + "metadata": {}, + "source": [ + "# Example workflow\n", + "\n", + "This notebook demonstrates downloading a short time series of GOES-R ABI imagery, merging the individual images into a single zarr dataset, and orthorectifying those images." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5c0d77eb-457b-4b4d-a669-73778a54637e", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/goes2go/data.py:665: FutureWarning: 'H' is deprecated and will be removed in a future version. Please use 'h' instead of 'H'.\n", + " within=pd.to_timedelta(config[\"nearesttime\"].get(\"within\", \"1h\")),\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/goes2go/NEW.py:188: FutureWarning: 'H' is deprecated and will be removed in a future version. Please use 'h' instead of 'H'.\n", + " within=pd.to_timedelta(config[\"nearesttime\"].get(\"within\", \"1h\")),\n" + ] + } + ], + "source": [ + "import goes_ortho as go\n", + "import xarray as xr\n", + "import geogif\n", + "import shutil" + ] + }, + { + "cell_type": "markdown", + "id": "0e01f57b-bbbe-4ebf-9c25-fff635326c2d", + "metadata": {}, + "source": [ + "First, specify the time range, location bounds, satellite, product (and if applicable, band and variable) that we'd like to access.\n", + "\n", + "We will also need to provide an API key for [OpenTopography.org](https://portal.opentopography.org/requestService?service=api) which you can create with a free account. This allows goes_ortho to access digital elevation models to perform the orthorectification step.\n", + "\n", + "The workflow below was developed to read a json file containing information about what we'd like to download. This was done to 1) allow these functions to run through github actions (still an experimental feature) and 2) keep a record of datasets we've downloaded. This is something that may change in the near future since it adds an unnecessary step for most use cases." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "28208c5f-e93c-4ce9-b4f7-c4c7ea14e12e", + "metadata": {}, + "outputs": [], + "source": [ + "# Make request file from user input\n", + "request_filepath = go.get_data.make_request_json(workflowName = \"example\",\n", + " startDatetime = \"2024-09-19T00:00:00Z\",\n", + " endDatetime = \"2024-09-20T00:59:00Z\",\n", + " bounds = go.get_data.bounds_from_geojson(\"grand_mesa.geojson\"),\n", + " satellite = \"goes18\",\n", + " product = \"ABI-L2-LSTC\",\n", + " band = 2,\n", + " variable = \"LST\",\n", + " apiKey = None,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "8eebc03e-2e78-44b8-beeb-e6914e454b2e", + "metadata": {}, + "source": [ + "The functions below demonstrate downloading GOES imagery using two different downloader packages: [goes2go](https://goes2go.readthedocs.io/en/latest/) and [goespy](https://github.com/spestana/goes-py) (the goespy functions are now integrated directly within the goes-ortho package). I have found goes2go is typically faster." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9fbfaf9c-938a-475a-a1f7-f2a36b062505", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Estimated 9 batches to download\n", + "Batch number 1\n", + "Download batch of imagery from 2024-09-19 00:00:00+00:00 to 2024-09-19 03:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:01<00:00, 1.81it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 2\n", + "Download batch of imagery from 2024-09-19 03:00:00+00:00 to 2024-09-19 06:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.92it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 3\n", + "Download batch of imagery from 2024-09-19 06:00:00+00:00 to 2024-09-19 09:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.86it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 4\n", + "Download batch of imagery from 2024-09-19 09:00:00+00:00 to 2024-09-19 12:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.92it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 5\n", + "Download batch of imagery from 2024-09-19 12:00:00+00:00 to 2024-09-19 15:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.88it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 6\n", + "Download batch of imagery from 2024-09-19 15:00:00+00:00 to 2024-09-19 18:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.85it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 7\n", + "Download batch of imagery from 2024-09-19 18:00:00+00:00 to 2024-09-19 21:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.78it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 8\n", + "Download batch of imagery from 2024-09-19 21:00:00+00:00 to 2024-09-20 00:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.78it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Batch number 9\n", + "Download batch of imagery from 2024-09-20 00:00:00+00:00 to 2024-09-20 03:00:00+00:00\n", + "📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].\n", + "Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 3/3 [00:00<00:00, 6.78it/s]" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Done\n", + "CPU times: user 5.34 s, sys: 465 ms, total: 5.8 s\n", + "Wall time: 13.5 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "%%time\n", + "filepaths = go.get_data.download_abi_goes2go(request_filepath)" + ] + }, + { + "cell_type": "markdown", + "id": "90be885d-c57f-4039-991c-9ef2543d3a78", + "metadata": {}, + "source": [ + "Merge all the individual NetCDFs that we just downloaded into a single zarr file." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "49523213-4a9a-4514-bed9-ef980387b572", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Dask dashboard at: http://127.0.0.1:8786/status\n", + "Workers: 6\n", + "Threads per worker: 2 \n", + "\n", + "read all the netcdf files, drop attributes where there are conflicts\n", + "configure chunking on variables with dimensions of (time, x, y)\n", + "Assign the dimensions of a chunk to variables to use for encoding afterwards\n", + "Assign the dimensions of a chunk to variables to use for encoding afterwards\n", + "Assign the dimensions of a chunk to variables to use for encoding afterwards\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/spestana/git/goes-ortho/src/goes_ortho/get_data.py:691: UserWarning: Times can't be serialized faithfully to int64 with requested units 'days since 2024-09-19T00:01:17'. Serializing with units 'hours since 2024-09-19T00:01:17' instead. Set encoding['dtype'] to floating point dtype to serialize with units 'days since 2024-09-19T00:01:17'. Set encoding['units'] to 'hours since 2024-09-19T00:01:17' to silence this warning .\n", + " ds.to_zarr(zarr_filepath)\n" + ] + } + ], + "source": [ + "zarr_filepath = f\"my_file.zarr\"\n", + "# remove if file already exists\n", + "shutil.rmtree(zarr_filepath, ignore_errors=True)\n", + "zarr_filepath = go.get_data.multi_nc_to_zarr(filepaths, zarr_filepath)" + ] + }, + { + "cell_type": "markdown", + "id": "226328c3-60f6-4126-84fa-1d4fe4dd2377", + "metadata": {}, + "source": [ + "Orthorectify the imagery in the zarr file" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "1c376545-ccc1-4774-b9ad-e22eeb7f55fe", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "https://portal.opentopography.org/API/globaldem?demtype=SRTMGL3&west=-108.368202&south=38.80429&east=-107.627676&north=39.211234&outputFormat=GTiff&API_Key=585b1d1639bc5ef8a4a5bdea7e45a8d1\n", + "/bin/gdalwarp -r cubic -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER -tr 30 30 -t_srs '+proj=lonlat +datum=GRS80' temp_SRTMGL3_DEM.tif temp_SRTMGL3_DEM_proj.tif\n", + "/bin/gdalwarp -r cubic -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER -tr 30 30 -t_srs '+proj=lonlat +datum=GRS80' temp_SRTMGL3_DEM.tif temp_SRTMGL3_DEM_proj.tif\n", + "Usage: gdalwarp [--help-general] [--formats]\n", + " [-s_srs srs_def] [-t_srs srs_def] [-to \"NAME=VALUE\"]* [-novshiftgrid]\n", + " [-order n | -tps | -rpc | -geoloc] [-et err_threshold]\n", + " [-refine_gcps tolerance [minimum_gcps]]\n", + " [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]\n", + " [-ovr level|AUTO|AUTO-n|NONE] [-wo \"NAME=VALUE\"] [-ot Byte/Int16/...] [-wt Byte/Int16]\n", + " [-srcnodata \"value [value...]\"] [-dstnodata \"value [value...]\"] -dstalpha\n", + " [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n", + " [-cutline datasource] [-cl layer] [-cwhere expression]\n", + " [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]\n", + " [-if format]* [-of format] [-co \"NAME=VALUE\"]* [-overwrite]\n", + " [-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]*\n", + " [-doo NAME=VALUE]*\n", + " srcfile* dstfile\n", + "\n", + "Available resampling methods:\n", + " near (default), bilinear, cubic, cubicspline, lanczos, average, rms,\n", + " mode, max, min, med, Q1, Q3, sum.\n", + "\n", + "RUNNING: make_ortho_map()\n", + "\n", + "Opening GOES ABI image...\n", + "\n", + "Get inputs: projection information from the ABI radiance product\n", + "...done\n", + "\n", + "Opening DEM file...\n", + "\n", + "Create 2D arrays of longitude and latitude from the DEM\n", + "...done\n", + "\n", + "For each grid cell in the DEM, compute the corresponding ABI scan angle (x and y, radians)\n", + "...done\n", + "\n", + "Create metadata dictionary about this map\n", + "...done\n", + "\n", + "Create pixel map dataset\n", + "...done\n", + "\n", + "Return the pixel map dataset.\n", + "\n", + "RUNNING: orthorectify_abi_rad()\n", + "\n", + "Does the projection info in the image match our mapping?\n", + "\n", + "Opening GOES ABI image...\t\t\tABI image value\tPixel map value\n", + "perspective_point_height + semi_major_axis:\t42164160.0\t42164160.0\n", + "semi_major_axis:\t\t\t\t6378137.0\t6378137.0\n", + "semi_minor_axis:\t\t\t\t6356752.31414\t6356752.31414\n", + "longitude_of_projection_origin:\t\t\t-137.0\t\t-137.0\n", + "...done\n", + "\n", + "Map (orthorectify) and clip the image to the pixel map for LST\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "ERROR 1: PROJ: proj_create: Error -9 (unknown elliptical parameter name)\n", + "ERROR 1: Translating source or target SRS failed:\n", + "+proj=lonlat +datum=GRS80\n", + "Child returned 1\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'netcdf4' fails while guessing\n", + " warnings.warn(f\"{engine!r} fails while guessing\", RuntimeWarning)\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'h5netcdf' fails while guessing\n", + " warnings.warn(f\"{engine!r} fails while guessing\", RuntimeWarning)\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'scipy' fails while guessing\n", + " warnings.warn(f\"{engine!r} fails while guessing\", RuntimeWarning)\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'netcdf4' fails while guessing\n", + " warnings.warn(f\"{engine!r} fails while guessing\", RuntimeWarning)\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'h5netcdf' fails while guessing\n", + " warnings.warn(f\"{engine!r} fails while guessing\", RuntimeWarning)\n", + "/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'scipy' fails while guessing\n", + " warnings.warn(f\"{engine!r} fails while guessing\", RuntimeWarning)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "...done\n", + "(27, 488, 889)\n", + "\n", + "Map (orthorectify) and clip the image to the pixel map for ABI Fixed Grid coordinates\n", + "...done\n", + "\n", + "Create zone labels for each unique pair of ABI Fixed Grid coordinates (for each orthorectified pixel footprint)\n", + "(488, 889)\n", + "...done\n", + "\n", + "Output this result to a new zarr file\n", + "Saving file as: test_ortho7.zarr\n", + "...done\n" + ] + } + ], + "source": [ + "bounds = go.get_data.bounds_from_geojson(\"grand_mesa.geojson\")\n", + "\n", + "api_key = \"585b1d1639bc5ef8a4a5bdea7e45a8d1\"\n", + "out_filename = 'test_ortho7.zarr'\n", + "shutil.rmtree(out_filename, ignore_errors=True)\n", + "go.orthorectify.ortho_zarr(\n", + " zarr_filepath,\n", + " ['LST'],\n", + " bounds,\n", + " api_key,\n", + " out_filename,\n", + " dem_filepath=None,\n", + " demtype=\"SRTMGL3\",\n", + " keep_dem=False,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "086756a9-9e42-4867-9410-10198ebdbe63", + "metadata": {}, + "source": [ + "Open the orthorectified imagery file" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "bb78ec87-2a3a-4418-a314-20ead2806a0d", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_zarr(out_filename)" + ] + }, + { + "cell_type": "markdown", + "id": "1ccb7d96-c2f8-40de-a626-3208fca8b0c0", + "metadata": {}, + "source": [ + "Make a gif animation to preview it" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "bd83d09e-a8a5-45b4-ba7d-7376e8db4cf4", + "metadata": {}, + "outputs": [], + "source": [ + "variable = 'LST'\n", + "# select our variable of interest\n", + "da = ds[variable]\n", + "\n", + "# create the gif animation\n", + "gif_bytes = geogif.dgif(\n", + " da,\n", + " fps=5,\n", + " cmap=\"Greys_r\",\n", + " date_format=\"%Y-%m-%d %H:%M:%S\",\n", + " date_position=\"ul\",\n", + " bytes=True,\n", + ").compute()\n", + "\n", + "# write gif to file\n", + "with open(f\"test6.gif\", \"wb\") as f:\n", + " f.write(gif_bytes)" + ] + }, + { + "cell_type": "markdown", + "id": "f15c4c93-9d86-4c88-a832-9f510ade26c2", + "metadata": {}, + "source": [ + "Take a look at the gif image we just made:\n", + "\n", + "\"GOES-18" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b9f5260-4d30-4782-8760-0a845484896d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "goes-test-env", + "language": "python", + "name": "goes-test-env" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/goes_ortho/get_data.py b/src/goes_ortho/get_data.py index 14a48d5..0c92e74 100644 --- a/src/goes_ortho/get_data.py +++ b/src/goes_ortho/get_data.py @@ -667,17 +667,26 @@ def multi_nc_to_zarr( verbose=True, ) as _: # read all the netcdf files, drop attributes where there are conflicts + print("read all the netcdf files, drop attributes where there are conflicts") ds = xr.open_mfdataset( - nc_filepaths, chunks={"time": 500}, combine_attrs="drop_conflicts" + nc_filepaths, + chunks={"time": 500}, + combine_attrs="drop_conflicts", + combine="by_coords", + compat="no_conflicts", ) # configure chunking on variables with dimensions of (time, x, y) + print("configure chunking on variables with dimensions of (time, x, y)") for variable, _ in ds.data_vars.items(): if ds[variable].dims == ("time", "y", "x"): ds[variable].data.rechunk( {0: -1, 1: "auto", 2: "auto"}, block_size_limit=1e8, balance=True ) # Assign the dimensions of a chunk to variables to use for encoding afterwards + print( + "Assign the dimensions of a chunk to variables to use for encoding afterwards" + ) t, y, x = ( ds[variable].data.chunks[0][0], ds[variable].data.chunks[1][0], diff --git a/src/goes_ortho/orthorectify.py b/src/goes_ortho/orthorectify.py index 51ffca7..a0dec4f 100644 --- a/src/goes_ortho/orthorectify.py +++ b/src/goes_ortho/orthorectify.py @@ -394,3 +394,350 @@ def ortho( os.remove(dem_filepath) return None + + +def make_ortho_map_zarr(goes_filepath, dem_filepath, out_filepath=None): + """ + For the entire DEM, determine the ABI scan angle coordinates for every DEM grid cell, taking into account the underlying terrain and satellite's viewing geometry. Create the mapping between GOES-R ABI pixels (zarr input file) and a DEM grid (geotiff input file) + + Parameters + ------------ + goes_filepath : str + filepath to GOES ABI zarr file + dem_filepath : str + filepath to digital elevation model (DEM), GeoTiff file + out_filepath : str + optional filepath and netcdf filename to save this map to, defaults to None + + Returns + ------------ + ds : xarray.Dataset + dataset of the map relating ABI Fixed Grid coordinates to latitude and longitude + + Examples + ------------ + + """ + + print("\nRUNNING: make_ortho_map()") + + # Open the GOES ABI image + print("\nOpening GOES ABI image...") + abi_image = xr.open_dataset(goes_filepath, decode_times=False) + # NOTE: for some reason (?) I sometimes get an error "ValueError: unable to decode time units 'seconds since 2000-01-01 12:00:00' with the default calendar. Try opening your dataset with decode_times=False." so I've added decode_times=False here. + # Get inputs: projection information from the ABI radiance product (values needed for geometry calculations) + print("\nGet inputs: projection information from the ABI radiance product") + req = abi_image.goes_imager_projection.semi_major_axis + rpol = abi_image.goes_imager_projection.semi_minor_axis + H = ( + abi_image.goes_imager_projection.perspective_point_height + + abi_image.goes_imager_projection.semi_major_axis + ) + lon_0 = abi_image.goes_imager_projection.longitude_of_projection_origin + e = 0.0818191910435 # GRS-80 eccentricity + print("...done") + + # Load DEM + print("\nOpening DEM file...") + dem = rioxarray.open_rasterio(dem_filepath) + dem = dem.where(dem != dem.attrs["_FillValue"])[0, :, :] # replace nodata with nans + dem = dem.fillna( + 0 + ) # fill nans with zeros for the ocean (temporary fix for fog project) + # dem = dem.where(dem!=0) # replace zeros with nans + # Create 2D arrays of longitude and latitude from the DEM + print("\nCreate 2D arrays of longitude and latitude from the DEM") + X, Y = np.meshgrid(dem.x, dem.y) # Lon and Lat of each DEM grid cell + Z = dem.values # elevation of each DEM grid cell + print("...done") + + # For each grid cell in the DEM, compute the corresponding ABI scan angle (x and y, radians) + print( + "\nFor each grid cell in the DEM, compute the corresponding ABI scan angle (x and y, radians)" + ) + abi_grid_x, abi_grid_y = LonLat2ABIangle(X, Y, Z, H, req, rpol, e, lon_0) + print("...done") + + # Create metadata dictionary about this map (should probably clean up metadata, adhere to some set of standards) + print("\nCreate metadata dictionary about this map") + metadata = { + # Information about the projection geometry: + "longitude_of_projection_origin": lon_0, + "semi_major_axis": req, + "semi_minor_axis": rpol, + "satellite_height": H, + "grs80_eccentricity": e, + "longitude_of_projection_origin_info": "longitude of geostationary satellite orbit", + "semi_major_axis_info": "semi-major axis of GRS 80 reference ellipsoid", + "semi_minor_axis_info": "semi-minor axis of GRS 80 reference ellipsoid", + "satellite_height_info": "distance from center of ellipsoid to satellite (perspective_point_height + semi_major_axis_info)", + "grs80_eccentricity_info": "eccentricity of GRS 80 reference ellipsoid", + # Information about the DEM source file + "dem_file": dem_filepath, + #'dem_crs' : dem.crs, + #'dem_transform' : dem.transform, + #'dem_res' : dem.res, + #'dem_ifov': -9999, # TO DO + "dem_file_info": "filename of dem file used to create this mapping", + "dem_crs_info": "coordinate reference system from DEM geotiff", + "dem_transform_info": "transform matrix from DEM geotiff", + "dem_res_info": "resolution of DEM geotiff", + "dem_ifov_info": "instantaneous field of view (angular size of DEM grid cell)", + # For each DEM grid cell, we have... + "dem_px_angle_x_info": "DEM grid cell X coordinate (east/west) scan angle in the ABI Fixed Grid", + "dem_px_angle_y_info": "DEM grid cell Y coordinate (north/south) scan angle in the ABI Fixed Grid", + "longitude_info": "longitude from DEM file", + "latitude_info": "latitude from DEM file", + "elevation_info": "elevation from DEM file", + } + print("...done") + + # Create pixel map dataset + print("\nCreate pixel map dataset") + ds = xr.Dataset( + {"elevation": (["latitude", "longitude"], dem.values)}, + coords={ + "longitude": (["longitude"], dem.x.data), + "latitude": (["latitude"], dem.y.data), + "dem_px_angle_x": (["latitude", "longitude"], abi_grid_x), + "dem_px_angle_y": (["latitude", "longitude"], abi_grid_y), + }, + attrs=metadata, + ) + # print(ds) + print("...done") + + if out_filepath is not None: + print("\nExport this pixel map along with the metadata (NetCDF with xarray)") + # Export this pixel map along with the metadata (NetCDF with xarray) + ds.to_netcdf(out_filepath, mode="w") + print("...done") + + # Return the pixel map dataset + print("\nReturn the pixel map dataset.") + + return ds + + +def orthorectify_abi_zarr(zarr_filepath, pixel_map, data_vars, out_filename=None): + """ + Using the pixel mapping for a specific ABI viewing geometry over a particular location, + orthorectify the ABI radiance values and return an xarray dataarray with those values. + + Parameters + ------------ + zarr_filepath : str + filepath to GOES ABI zarr file + pixel_map : xarray.Dataset + dataset of the map relating ABI Fixed Grid coordinates to latitude and longitude + data_vars : list + list of variable names from the GOES ABI NetCDF file we wish to extract + out_filename : str + optional filepath and zarr filename to save the orthorectified image to, defaults to None + + Returns + ------------ + pixel_map : xarray.Dataset + dataset of the orthorectified GOES ABI image + + Examples + ------------ + + """ + print("\nRUNNING: orthorectify_abi_rad()") + + # First check, Does the projection info in the image match our mapping? + print("\nDoes the projection info in the image match our mapping?") + # Open the GOES ABI image + print("\nOpening GOES ABI image...\t\t\tABI image value\tPixel map value") + abi_image = xr.open_dataset(zarr_filepath) + + print( + "perspective_point_height + semi_major_axis:\t{}\t{}".format( + abi_image.goes_imager_projection.perspective_point_height + + abi_image.goes_imager_projection.semi_major_axis, + pixel_map.satellite_height, + ) + ) + print( + "semi_major_axis:\t\t\t\t{}\t{}".format( + abi_image.goes_imager_projection.semi_major_axis, pixel_map.semi_major_axis + ) + ) + print( + "semi_minor_axis:\t\t\t\t{}\t{}".format( + abi_image.goes_imager_projection.semi_minor_axis, pixel_map.semi_minor_axis + ) + ) + print( + "longitude_of_projection_origin:\t\t\t{}\t\t{}".format( + abi_image.goes_imager_projection.longitude_of_projection_origin, + pixel_map.longitude_of_projection_origin, + ) + ) + print("...done") + + # Map (orthorectify) and clip the image to the pixel map for each data variable we want + for var in data_vars: + print( + "\nMap (orthorectify) and clip the image to the pixel map for {}".format( + var + ) + ) + abi_var_values = abi_image.sel( + x=pixel_map.dem_px_angle_x, y=pixel_map.dem_px_angle_y, method="nearest" + )[var].values + print("...done") + + # Create a new xarray dataset with the orthorectified ABI radiance values, + # Lat, Lon, Elevation, and metadata from the pixel map. + print(abi_var_values.shape) + pixel_map[var] = (["time", "latitude", "longitude"], abi_var_values) + + ## If we are looking at an ABI-L1b-Rad product, create either a reflectance (bands 1-6) or brightness temperautre (bands 7-16) dataset + # if var == "Rad": + # # if we are looking at bands 1-6, compute reflectance + # if abi_image.band_id.values[0] <= 6: + # pixel_map["ref"] = goesReflectance( + # pixel_map[var], abi_image.kappa0.values + # ) + # # else, compute brightness temperature for bands 7-16 + # else: + # pixel_map["tb"] = goesBrightnessTemp( + # pixel_map[var], + # abi_image.planck_fk1.values, + # abi_image.planck_fk2.values, + # abi_image.planck_bc1.values, + # abi_image.planck_bc2.values, + # ) + + # Map (orthorectify) the original ABI Fixed Grid coordinate values to the new pixels for reference + print( + "\nMap (orthorectify) and clip the image to the pixel map for ABI Fixed Grid coordinates" + ) + abi_fixed_grid_x_values = abi_image.sel( + x=pixel_map.dem_px_angle_x.values.ravel(), method="nearest" + ).x.values + abi_fixed_grid_y_values = abi_image.sel( + y=pixel_map.dem_px_angle_y.values.ravel(), method="nearest" + ).y.values + abi_fixed_grid_x_values_reshaped = np.reshape( + abi_fixed_grid_x_values, pixel_map.dem_px_angle_x.shape + ) + abi_fixed_grid_y_values_reshaped = np.reshape( + abi_fixed_grid_y_values, pixel_map.dem_px_angle_y.shape + ) + pixel_map["abi_fixed_grid_x"] = ( + ("latitude", "longitude"), + abi_fixed_grid_x_values_reshaped, + ) + pixel_map["abi_fixed_grid_y"] = ( + ("latitude", "longitude"), + abi_fixed_grid_y_values_reshaped, + ) + print("...done") + + # drop DEM from dataset + # pixel_map = pixel_map.drop(['elevation']) + + print( + "\nCreate zone labels for each unique pair of ABI Fixed Grid coordinates (for each orthorectified pixel footprint)" + ) + # Found this clever solution here: https://stackoverflow.com/a/32326297/11699349 + # Create unique values for every "zone" (the GOES ABI pixel footprints) with the same ABI Fixed Grid X and Y values + unique_values = ( + pixel_map.abi_fixed_grid_x.values + * (pixel_map.abi_fixed_grid_y.values.max() + 1) + + pixel_map.abi_fixed_grid_y.values + ) + # Find the index of all unique values we just created + _, idx = np.unique(unique_values, return_inverse=True) + # Use these indices, reshaped to the original shape, as our zone labels + zone_labels = idx.reshape(pixel_map.abi_fixed_grid_y.values.shape) + print(zone_labels.shape) + # Add the zone_labels to the dataset + pixel_map["zone_labels"] = (("latitude", "longitude"), zone_labels) + print("...done") + + # add back in the time coordinates + pixel_map = pixel_map.assign_coords({"time": abi_image.time.values}) + pixel_map = pixel_map.reset_coords(drop=True) + + # Output this result to a new NetCDF file + print("\nOutput this result to a new zarr file") + if out_filename is None: + out_filename = abi_image.dataset_name + "_ortho.zarr" + print("Saving file as: {}".format(out_filename)) + + pixel_map.to_zarr(out_filename) + print("...done") + + return pixel_map + + +def ortho_zarr( + goes_image_path, + data_vars, + bounds, + api_key, + new_goes_filename, + dem_filepath=None, + demtype="SRTMGL3", + keep_dem=True, +): + """ + Wraps around get_dem(), make_ortho_map_zarr(), orthorectify_abi_zarr() + + Parameters + ------------ + goes_image_path : str + filepath to GOES ABI zarr file + data_vars : list + list of variable names from the GOES ABI zarr file we wish to extract + bounds : list + longitude and latitude bounds to clip and orthorectify GOES ABI image, like [min_lon, min_lat, max_lon, max_lat] + api_key : str + Opentopography.org API key, can be created at https://portal.opentopography.org/requestService?service=api + new_goes_filename : str + new filepath and filename to save the orthorectified image to + dem_filepath : str + filepath to save DEM to, defaults to None + demtype : str + DEM from Opentopography.org, see documentation in get_data.get_dem() + keep_dem : bool + option to save DEM file or delete after use + + Returns + ------------ + None + + Examples + ------------ + + """ + + if dem_filepath is None: + dem_filepath = "temp_{demtype}_DEM.tif".format( + demtype=demtype, + ) + get_dem( + demtype=demtype, + bounds=bounds, + api_key=api_key, + out_fn=dem_filepath, + proj="+proj=lonlat +datum=GRS80", + ) # make sure to convert to GRS80 ellipsoid model GOES ABI fixed grid uses + + # create the mapping between scan angle coordinates and lat/lon given the GOES satellite position and our DEM + goes_ortho_map = make_ortho_map_zarr(goes_image_path, dem_filepath) + + # Apply the "ortho map" and save a new NetCDF file with data variables from the original file + _ = orthorectify_abi_zarr( + goes_image_path, goes_ortho_map, data_vars, out_filename=new_goes_filename + ) + + # If keep_dem is False, delete the temporary DEM file we downloaded + if not keep_dem: + os.remove(dem_filepath) + + return None From 405788f6067dab0b299063c07148060a9d8eb15c Mon Sep 17 00:00:00 2001 From: spestana Date: Fri, 27 Sep 2024 02:34:06 +0000 Subject: [PATCH 3/3] Update toc with new examples --- ...-workflow.ipynb => example_workflow.ipynb} | 0 docs/index.rst | 19 +++++++++++-------- 2 files changed, 11 insertions(+), 8 deletions(-) rename docs/examples/{example-workflow.ipynb => example_workflow.ipynb} (100%) diff --git a/docs/examples/example-workflow.ipynb b/docs/examples/example_workflow.ipynb similarity index 100% rename from docs/examples/example-workflow.ipynb rename to docs/examples/example_workflow.ipynb diff --git a/docs/index.rst b/docs/index.rst index b0f5e7e..88604c2 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -8,14 +8,7 @@ installation examples/quick_start examples/download - examples/subset_abi_netcdf_example - examples/orthorectify_abi_example - examples/orthorectify_abi_example2 - examples/ortho-horizontal-offset - examples/make_abi_timeseries_example - examples/jsontest - examples/goes-orthorectify - + examples/example_workflow .. toctree:: :caption: API Reference @@ -27,3 +20,13 @@ geometry clip rad + +.. toctree:: + :caption: Old Examples + :maxdepth: 1 + examples/orthorectify_abi_example + examples/orthorectify_abi_example2 + examples/ortho-horizontal-offset + examples/make_abi_timeseries_example + examples/jsontest + examples/goes-orthorectify