CLAAS Grid Inconsistency

claas_grid_inconsistency.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shift_area(area, shift=1500):\n",
    "    \"\"\"Shift area extent by roughly half a pixel to S-E to account for the image offset.\"\"\"\n",
    "    llx, lly, urx, ury = area.area_extent\n",
    "    return area.copy(area_extent=[llx+shift, lly-shift, urx+shift, ury-shift])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_true_claas_extent(full_disk_area_def):\n",
    "    \"\"\"Get true area extent of CLAAS grid by subsetting the full disk.\"\"\"\n",
    "    full_disk_size = 3712\n",
    "    claas_size = 3636\n",
    "    offset = (full_disk_size - claas_size + 1) // 2\n",
    "    x, y = full_disk_area_def.get_proj_vectors()\n",
    "    x_claas = x[offset:-offset]\n",
    "    y_claas = y[::-1][offset:-offset]\n",
    "    return np.array([x_claas[0], y_claas[0], x_claas[-1], y_claas[-1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "...lib/python3.7/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems\n",
      "  proj_string = self.to_proj4()\n"
     ]
    }
   ],
   "source": [
    "from pyresample.geometry import AreaDefinition\n",
    "\n",
    "# Define full disk area definitions.\n",
    "\n",
    "MAJOR_AXIS_OF_EARTH_ELLIPSOID = 6378169.0\n",
    "MINOR_AXIS_OF_EARTH_ELLIPSOID = 6356583.8\n",
    "SATELLITE_ALTITUDE = 35785831.0\n",
    "PROJECTION_LONGITUDE = 0.0\n",
    "PROJ_DICT = {\n",
    "    \"a\": MAJOR_AXIS_OF_EARTH_ELLIPSOID,\n",
    "    \"b\": MINOR_AXIS_OF_EARTH_ELLIPSOID,\n",
    "    \"h\": SATELLITE_ALTITUDE,\n",
    "    \"lon_0\": PROJECTION_LONGITUDE,\n",
    "    \"proj\": \"geos\",\n",
    "    \"units\": \"m\",\n",
    "}\n",
    "\n",
    "# 1. No offset. This is identical to \"msg_seviri_fes_3km\" in satpy/etc/areas.yaml\n",
    "full_disk_no_offset = AreaDefinition(\n",
    "    area_id=\"full_disk_no_offset\",\n",
    "    description=\"full_disk_no_offset\",\n",
    "    proj_id=\"full_disk_no_offset\",\n",
    "    projection=PROJ_DICT,\n",
    "    area_extent=[-5570248.686685662, -5567248.28340708, 5567248.28340708, 5570248.686685662],\n",
    "    width=3712,\n",
    "    height=3712,\n",
    ")\n",
    "\n",
    "# 2. With offset\n",
    "full_disk_with_offset = shift_area(full_disk_no_offset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Area extents from CLAAS netcdf files\n",
    "extent_claas2_nc = np.array([-5456233.41938636, -5453233.01608472, 5453233.01608472, 5456233.41938636])\n",
    "extent_claas21_nc = np.array([-5454733.16046029, -5451732.75718171, 5451732.75718171, 5454733.16046029])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1500.25892607, -1500.25890301,  1500.25890301,  1500.25892607])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CLAAS-2 area extent is larger by one pixel... Shoot!\n",
    "extent_claas2_nc - extent_claas21_nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1500.25892607 -1500.25890301  1500.25890301  1500.25892607]\n",
      "[-3.00025893e+03 -2.58903010e-01  2.58903010e-01  3.00025893e+03]\n"
     ]
    }
   ],
   "source": [
    "# Differences between CLAAS-2 grid and true grids\n",
    "print(extent_claas2_nc - get_true_claas_extent(full_disk_no_offset))\n",
    "print(extent_claas2_nc - get_true_claas_extent(full_disk_with_offset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.79396772e-09  9.31322575e-10 -9.31322575e-10 -2.79396772e-09]\n",
      "[-1500.  1500. -1500.  1500.]\n"
     ]
    }
   ],
   "source": [
    "# Differences between CLAAS-2.1 grid and true grids\n",
    "print(extent_claas21_nc - get_true_claas_extent(full_disk_no_offset))\n",
    "print(extent_claas21_nc - get_true_claas_extent(full_disk_with_offset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "girafe",
   "language": "python",
   "name": "girafe"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}