Example workflow#
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.
[1]:
import goes_ortho as go
import xarray as xr
import geogif
import shutil
/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'.
within=pd.to_timedelta(config["nearesttime"].get("within", "1h")),
/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'.
within=pd.to_timedelta(config["nearesttime"].get("within", "1h")),
First, specify the time range, location bounds, satellite, product (and if applicable, band and variable) that we’d like to access.
We will also need to provide an API key for OpenTopography.org which you can create with a free account. This allows goes_ortho to access digital elevation models to perform the orthorectification step.
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.
[2]:
# Make request file from user input
request_filepath = go.get_data.make_request_json(workflowName = "example",
startDatetime = "2024-09-19T00:00:00Z",
endDatetime = "2024-09-20T00:59:00Z",
bounds = go.get_data.bounds_from_geojson("grand_mesa.geojson"),
satellite = "goes18",
product = "ABI-L2-LSTC",
band = 2,
variable = "LST",
apiKey = None,
)
The functions below demonstrate downloading GOES imagery using two different downloader packages: goes2go and goespy (the goespy functions are now integrated directly within the goes-ortho package). I have found goes2go is typically faster.
[3]:
%%time
filepaths = go.get_data.download_abi_goes2go(request_filepath)
Estimated 9 batches to download
Batch number 1
Download batch of imagery from 2024-09-19 00:00:00+00:00 to 2024-09-19 03:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:01<00:00, 1.81it/s]
Batch number 2
Download batch of imagery from 2024-09-19 03:00:00+00:00 to 2024-09-19 06:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.92it/s]
Batch number 3
Download batch of imagery from 2024-09-19 06:00:00+00:00 to 2024-09-19 09:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.86it/s]
Batch number 4
Download batch of imagery from 2024-09-19 09:00:00+00:00 to 2024-09-19 12:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.92it/s]
Batch number 5
Download batch of imagery from 2024-09-19 12:00:00+00:00 to 2024-09-19 15:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.88it/s]
Batch number 6
Download batch of imagery from 2024-09-19 15:00:00+00:00 to 2024-09-19 18:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.85it/s]
Batch number 7
Download batch of imagery from 2024-09-19 18:00:00+00:00 to 2024-09-19 21:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.78it/s]
Batch number 8
Download batch of imagery from 2024-09-19 21:00:00+00:00 to 2024-09-20 00:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.78it/s]
Batch number 9
Download batch of imagery from 2024-09-20 00:00:00+00:00 to 2024-09-20 03:00:00+00:00
📦 Finished downloading [3] files to [/home/spestana/data/noaa-goes18/ABI-L2-LSTC].
Cropping image batch to [-108.368202, 38.80429, -107.627676, 39.211234]
100%|██████████| 3/3 [00:00<00:00, 6.78it/s]
Done
CPU times: user 5.34 s, sys: 465 ms, total: 5.8 s
Wall time: 13.5 s
Merge all the individual NetCDFs that we just downloaded into a single zarr file.
[7]:
zarr_filepath = f"my_file.zarr"
# remove if file already exists
shutil.rmtree(zarr_filepath, ignore_errors=True)
zarr_filepath = go.get_data.multi_nc_to_zarr(filepaths, zarr_filepath)
Dask dashboard at: http://127.0.0.1:8786/status
Workers: 6
Threads per worker: 2
read all the netcdf files, drop attributes where there are conflicts
configure chunking on variables with dimensions of (time, x, y)
Assign the dimensions of a chunk to variables to use for encoding afterwards
Assign the dimensions of a chunk to variables to use for encoding afterwards
Assign the dimensions of a chunk to variables to use for encoding afterwards
/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 .
ds.to_zarr(zarr_filepath)
Orthorectify the imagery in the zarr file
[18]:
bounds = go.get_data.bounds_from_geojson("grand_mesa.geojson")
api_key = "585b1d1639bc5ef8a4a5bdea7e45a8d1"
out_filename = 'test_ortho7.zarr'
shutil.rmtree(out_filename, ignore_errors=True)
go.orthorectify.ortho_zarr(
zarr_filepath,
['LST'],
bounds,
api_key,
out_filename,
dem_filepath=None,
demtype="SRTMGL3",
keep_dem=False,
)
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
/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
/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
Usage: gdalwarp [--help-general] [--formats]
[-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]* [-novshiftgrid]
[-order n | -tps | -rpc | -geoloc] [-et err_threshold]
[-refine_gcps tolerance [minimum_gcps]]
[-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
[-ovr level|AUTO|AUTO-n|NONE] [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
[-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha
[-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
[-cutline datasource] [-cl layer] [-cwhere expression]
[-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
[-if format]* [-of format] [-co "NAME=VALUE"]* [-overwrite]
[-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]*
[-doo NAME=VALUE]*
srcfile* dstfile
Available resampling methods:
near (default), bilinear, cubic, cubicspline, lanczos, average, rms,
mode, max, min, med, Q1, Q3, sum.
RUNNING: make_ortho_map()
Opening GOES ABI image...
Get inputs: projection information from the ABI radiance product
...done
Opening DEM file...
Create 2D arrays of longitude and latitude from the DEM
...done
For each grid cell in the DEM, compute the corresponding ABI scan angle (x and y, radians)
...done
Create metadata dictionary about this map
...done
Create pixel map dataset
...done
Return the pixel map dataset.
RUNNING: orthorectify_abi_rad()
Does the projection info in the image match our mapping?
Opening GOES ABI image... ABI image value Pixel map value
perspective_point_height + semi_major_axis: 42164160.0 42164160.0
semi_major_axis: 6378137.0 6378137.0
semi_minor_axis: 6356752.31414 6356752.31414
longitude_of_projection_origin: -137.0 -137.0
...done
Map (orthorectify) and clip the image to the pixel map for LST
ERROR 1: PROJ: proj_create: Error -9 (unknown elliptical parameter name)
ERROR 1: Translating source or target SRS failed:
+proj=lonlat +datum=GRS80
Child returned 1
/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'netcdf4' fails while guessing
warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'h5netcdf' fails while guessing
warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'scipy' fails while guessing
warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'netcdf4' fails while guessing
warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'h5netcdf' fails while guessing
warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/home/spestana/.conda/envs/goes-test-env/lib/python3.12/site-packages/xarray/backends/plugins.py:149: RuntimeWarning: 'scipy' fails while guessing
warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
...done
(27, 488, 889)
Map (orthorectify) and clip the image to the pixel map for ABI Fixed Grid coordinates
...done
Create zone labels for each unique pair of ABI Fixed Grid coordinates (for each orthorectified pixel footprint)
(488, 889)
...done
Output this result to a new zarr file
Saving file as: test_ortho7.zarr
...done
Open the orthorectified imagery file
[19]:
ds = xr.open_zarr(out_filename)
Make a gif animation to preview it
[22]:
variable = 'LST'
# select our variable of interest
da = ds[variable]
# create the gif animation
gif_bytes = geogif.dgif(
da,
fps=5,
cmap="Greys_r",
date_format="%Y-%m-%d %H:%M:%S",
date_position="ul",
bytes=True,
).compute()
# write gif to file
with open(f"test6.gif", "wb") as f:
f.write(gif_bytes)
Take a look at the gif image we just made:
[ ]: