Reprojecting GOES16 data for conversion to COG with rioxarray

goes16-rioxarray-reprojection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a014e2b-fee7-4670-8ae2-715c03fde690",
   "metadata": {},
   "source": [
    "# Load and reproject GOES16 data with rioxarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "68586e39-7cb0-4803-8cbe-9069350aa2cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import s3fs \n",
    "import xarray as xr\n",
    "import rioxarray\n",
    "from pyproj import CRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "10624d5f-5de4-4dce-96fe-4290768c1570",
   "metadata": {},
   "outputs": [],
   "source": [
    "def goes_to_wgs(ds, variable_name):\n",
    "\n",
    "    # convert coordinates to meters https://gis.stackexchange.com/a/350006\n",
    "    sat_height = ds[\"goes_imager_projection\"].attrs[\"perspective_point_height\"]\n",
    "    ds = ds.assign_coords({\n",
    "        \"x\": ds[\"x\"].values * sat_height,\n",
    "        \"y\": ds[\"y\"].values * sat_height,\n",
    "    })\n",
    "    \n",
    "    crs = CRS.from_cf(ds[\"goes_imager_projection\"].attrs)\n",
    "    ds.rio.write_crs(crs.to_string(), inplace=True)\n",
    "    \n",
    "    da = ds[variable_name]\n",
    "\n",
    "    return da.rio.reproject(\"epsg:4326\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f90c18c-f087-4bd8-b2b3-0ded44696b82",
   "metadata": {},
   "source": [
    "## Straight from AWS S3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "20167f1c-b5c1-4efc-8240-3b5443205904",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = s3fs.S3FileSystem(anon=True)\n",
    "\n",
    "path = \"s3://noaa-goes16/ABI-L1b-RadF/2017/129/11/OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc\"\n",
    "\n",
    "with fs.open(path) as fileObj:\n",
    "    with xr.open_dataset(fileObj, engine='h5netcdf') as ds:\n",
    "\n",
    "        da = goes_to_wgs(ds, variable_name=\"Rad\")\n",
    "    \n",
    "        COG_PROFILE = {\"driver\": \"COG\", \"compress\": \"DEFLATE\", \"predictor\": 2}\n",
    "        da.rio.to_raster(\"goes16.tif\", **COG_PROFILE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baf3831d-f12c-4df5-ae3b-667569b054fa",
   "metadata": {},
   "source": [
    "## Local file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e7762562-9fe8-4355-9eff-0a87114d1095",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 50.2M  100 50.2M    0     0  11.7M      0  0:00:04  0:00:04 --:--:-- 11.7M0:00:03 9166k\n"
     ]
    }
   ],
   "source": [
    "!curl -O https://noaa-goes16.s3.amazonaws.com/ABI-L1b-RadF/2017/129/11/OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "edceea29-16e4-45cc-966b-3c7e3969d8e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = \"OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc\"\n",
    "\n",
    "with xr.open_dataset(path, engine='h5netcdf') as ds:\n",
    "\n",
    "    da = goes_to_wgs(ds, variable_name=\"Rad\")\n",
    "\n",
    "    COG_PROFILE = {\"driver\": \"COG\", \"compress\": \"DEFLATE\", \"predictor\": 2}\n",
    "    da.rio.to_raster(\"goes16.tif\", **COG_PROFILE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce9b32fe-59d9-46fd-8b28-961943cb2758",
   "metadata": {},
   "source": [
    "## With GDAL CLI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b92dcd90-f1bf-4b88-9cfc-db3bfa6f7095",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning 1: Unhandled X/Y axis unit rad. SRS will ignore axis unit and be likely wrong.\n",
      "Input file size is 10848, 10848\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "!gdal_translate -ot float32 -unscale -CO COMPRESS=deflate NETCDF:\"{path}\":Rad fulldisk.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "49273e91-f1c2-473c-aac7-f26d083363e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating output file that is 11178P x 10507L.\n",
      "Using internal nodata values (e.g. 1023) for image fulldisk.tif.\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "!gdalwarp -t_srs EPSG:4326 -dstnodata -999.0 -of COG fulldisk.tif fulldisk_geo.tif"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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
}