03_match_projection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Match map projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Requirements:\n",
    "> conda install -c conda-forge gdal=3 jupyter matplotlib rasterio tuiview\n",
    "\n",
    "**Data for lesson**\n",
    "\n",
    "CTX Image: \n",
    "\n",
    "https://astrogeology.usgs.gov/search/map/Mars/Mars2020/landing_site/J03_045994_1986_J03_046060_1986_20m_DTM\n",
    "\n",
    "This file is in GeoTIFF format with optional ISIS3 and PDS3 labels available for download. Copy into a \"ctx/\" directory below this notebook.\n",
    "\n",
    "files:\n",
    "\n",
    "* J03_045994_1986_J03_046060_1986_20m_DTM.tif\n",
    "\n",
    "---\n",
    "\n",
    "HiRISE Images:\n",
    "https://www.uahirise.org/dtm/dtm.php?ID=ESP_045994_1985\n",
    "\n",
    "This DEM/DTM files is in PDS3 format. Copy into a \"hirise/\" directory below this notebook.\n",
    "\n",
    "files:\n",
    "\n",
    "* DTEEC_045994_1985_046060_1985_U01.IMG\n",
    "\n",
    "---\n",
    "\n",
    "Note: \n",
    "* Using **!** allows one to run a routine outside of the notebook in the terminal\n",
    "* $variable, sends the variable's value to the terminal (command-line)\n",
    "* \"//\" for path works in Linux (Mac) and Windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs=True\n"
     ]
    }
   ],
   "source": [
    "match_to = \"ctx/J03_045994_1986_J03_046060_1986_20m_DTM.tif\"\n",
    "\n",
    "import rasterio\n",
    "t_srs = rasterio.open(match_to).crs.to_proj4()\n",
    "print(t_srs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=eqc +lat_ts=15 +lat_0=0 +lon_0=77.44 +x_0=0 +y_0=0 +R=3394839.8133163 +units=m +no_defs=True\n",
      "\n",
      "Creating output file that is 6864P x 14175L.\n",
      "Processing hirise//DTEEC_045994_1985_046060_1985_U01.IMG [1/1] : 0Using internal nodata values (e.g. -3.40282e+38) for image hirise//DTEEC_045994_1985_046060_1985_U01.IMG.\n",
      "Copying nodata values from source hirise//DTEEC_045994_1985_046060_1985_U01.IMG to destination hirise//DTEEC_045994_1985_046060_1985_U01.IMGout_srs.tif.\n",
      "...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "replica = \"hirise//DTEEC_045994_1985_046060_1985_U01.IMG\" \n",
    "to = replica + \"out_srs.tif\"\n",
    "\n",
    "# what are we map projecting from?\n",
    "r_srs = rasterio.open(replica).crs.to_proj4()\n",
    "print(r_srs + \"\\n\")\n",
    "#note the localized Mars Radius, common in HiRISE labels\n",
    "\n",
    "!gdalwarp -overwrite -t_srs \"$t_srs\" $replica $to -r bilinear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs=True\n"
     ]
    }
   ],
   "source": [
    "# check the output projection matches the \"match_to\" projection\n",
    "\n",
    "t_srs = rasterio.open(to).crs.to_proj4()\n",
    "print(t_srs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a oneliner if used in bash:\n",
    "> gdalwarp -t_srs < ( gdalsrsinfo -o wkt match_to.tif ) replica.tif  to.tif\n",
    "\n",
    "from: https://gis.stackexchange.com/questions/239013/reprojecting-a-raster-to-match-another-raster-using-gdal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here tuiview will prove these teo images are not in the same space.\n",
    "# tuiview does not support projection-on-fly (like QGIS and ArcGIS)\n",
    "#!tuiview $match_to $replica >/dev/null 2>&1\n",
    "\n",
    "# Here tuiview reports the projects are different, but they are the same\n",
    "# and thus the images are now co-registered in tuiview\n",
    "!tuiview $match_to $to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Let's try a short Python script to match the projection and **force the image size (and thus the output resolution)** to be the same. See match_projection.py code list on bottom of this notebook.\n",
    "\n",
    "from: https://gis.stackexchange.com/questions/234512/matching-two-rasters-with-different-projections-and-resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from match_projection import match_projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nodata: -3.4028226550889045e+38\n",
      "band:  1\n",
      "output file create: hirise//DTEEC_045994_1985_046060_1985_U01.IMG_match.tif\n"
     ]
    }
   ],
   "source": [
    "match_to = \"ctx//J03_045994_1986_J03_046060_1986_20m_DTM.tif\"\n",
    "replica = \"hirise//DTEEC_045994_1985_046060_1985_U01.IMG\" \n",
    "\n",
    "out = match_projection(match_to, replica)\n",
    "print(\"output file create: \" + out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs=True\n"
     ]
    }
   ],
   "source": [
    "o_srs = rasterio.open(to).crs.to_proj4()\n",
    "print(o_srs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tuiview $match_to $out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus - GDAL difference\n",
    "\n",
    "gdal_calc.py - https://gdal.org/programs/gdal_calc.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# limitation - gdal_calc.py must have files of the same size and \n",
    "#   it also doesn't check projections\n",
    "\n",
    "diff = \"hirise//out_diff.tif\"\n",
    "\n",
    "#!gdal_calc.py\n",
    "!gdal_calc.py --quiet --calc \"A - B\" --outfile $diff -A $match_to -B $out \n",
    "print(\"created: \" + diff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tuiview $diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### match_projection.py code listing\n",
    "   \n",
    "original code (updated) from: https://github.com/jgomezdans/eoldas_ng_observations/blob/master/eoldas_ng_observations/eoldas_observation_helpers.py#L29\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from osgeo import gdal\n",
    "\n",
    "def match_projection ( primary, replica ):\n",
    "    \"\"\"This function reprojects an image (``replica``) to\n",
    "    match the extent, resolution and projection of another\n",
    "    (``primary``) using GDAL. The newly reprojected image\n",
    "    is a GDAL GTIFF file.\n",
    "\n",
    "    Parameters\n",
    "    -------------\n",
    "    primary: str \n",
    "        A filename (with full path if required) with the \n",
    "        primary image (that that will be taken as a reference\n",
    "        and matched)\n",
    "    replica: str \n",
    "        A filename (with path if needed) with the image\n",
    "        that will be reprojected\n",
    "    Returns\n",
    "    ----------\n",
    "    The reprojected filename\n",
    "    \"\"\"\n",
    "    replica_ds = gdal.Open( replica )\n",
    "    if (replica_ds is None):\n",
    "       raise (IOError, \"GDAL could not open replica file %s \" \\\n",
    "            % replica)\n",
    "    replica_proj = replica_ds.GetProjection()\n",
    "    replica_geotrans = replica_ds.GetGeoTransform()\n",
    "    data_type = replica_ds.GetRasterBand(1).DataType\n",
    "    data_nodata = replica_ds.GetRasterBand(1).GetNoDataValue()\n",
    "    n_bands = replica_ds.RasterCount\n",
    "\n",
    "    primary_ds = gdal.Open( primary )\n",
    "    if (primary_ds is None):\n",
    "        raise (IOError, \"GDAL could not open primary file %s \" \\\n",
    "            % primary)\n",
    "    primary_proj = primary_ds.GetProjection()\n",
    "    primary_geotrans = primary_ds.GetGeoTransform()\n",
    "    w = primary_ds.RasterXSize\n",
    "    h = primary_ds.RasterYSize\n",
    "\n",
    "    dst_filename = replica + \"_match.tif\"\n",
    "    dst_ds = gdal.GetDriverByName(\"GTiff\").Create(dst_filename,\n",
    "                                                w, h, n_bands, data_type)\n",
    "    dst_ds.SetGeoTransform( primary_geotrans )\n",
    "    dst_ds.SetProjection( primary_proj)\n",
    "    print (\"nodata: \" + str(data_nodata))\n",
    "    \n",
    "    for i in range(1, n_bands + 1):\n",
    "        if data_nodata is not None:\n",
    "            dst_ds.GetRasterBand(i).SetNoDataValue(data_nodata)\n",
    "            print (\"band:  \" + str(i))\n",
    "\n",
    "    gdal.ReprojectImage( replica_ds, dst_ds, replica_proj,\n",
    "                         primary_proj, gdal.GRA_Bilinear)\n",
    "    dst_ds = None  # Flush to disk\n",
    "    return (dst_filename)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}