ravin_dem_difference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from osgeo import gdal, gdalconst\n",
    "import numpy as np\n",
    "\n",
    "def align_dems(base_dem_path, target_dem_path, aligned_dem_path):\n",
    "    # Open base DEM (reference) and get its properties\n",
    "    base_dem = gdal.Open(base_dem_path, gdalconst.GA_ReadOnly)\n",
    "    base_dem_proj = base_dem.GetProjection()\n",
    "    base_dem_geotrans = base_dem.GetGeoTransform()\n",
    "    base_dem_xsize = base_dem.RasterXSize\n",
    "    base_dem_ysize = base_dem.RasterYSize\n",
    "    base_dem_nodata = base_dem.GetRasterBand(1).GetNoDataValue()\n",
    "\n",
    "    # Open the target DEM\n",
    "    target_dem = gdal.Open(target_dem_path, gdalconst.GA_ReadOnly)\n",
    "    target_dem_nodata = target_dem.GetRasterBand(1).GetNoDataValue()\n",
    "\n",
    "    # Create a new aligned DEM file\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    aligned_dem = driver.Create(aligned_dem_path, base_dem_xsize, base_dem_ysize, 1, target_dem.GetRasterBand(1).DataType)\n",
    "    aligned_dem.SetProjection(base_dem_proj)\n",
    "    aligned_dem.SetGeoTransform(base_dem_geotrans)\n",
    "\n",
    "    # Set NoData value for aligned DEM\n",
    "    nodata_value = target_dem_nodata if target_dem_nodata is not None else base_dem_nodata\n",
    "    aligned_dem.GetRasterBand(1).SetNoDataValue(nodata_value)\n",
    "\n",
    "    # Perform the alignment/resampling\n",
    "    gdal.ReprojectImage(target_dem, aligned_dem, target_dem.GetProjection(), base_dem_proj, gdalconst.GRA_Bilinear)\n",
    "\n",
    "    # Close the datasets\n",
    "    base_dem = None\n",
    "    target_dem = None\n",
    "    aligned_dem = None\n",
    "\n",
    "    # Modify the aligned DEM to set 0 values to NoData\n",
    "    modify_dem(aligned_dem_path, nodata_value)\n",
    "\n",
    "def modify_dem(dem_path, nodata_value):\n",
    "    # Open the aligned DEM\n",
    "    dem = gdal.Open(dem_path, gdalconst.GA_Update)\n",
    "    band = dem.GetRasterBand(1)\n",
    "    data = band.ReadAsArray()\n",
    "\n",
    "    # Set 0 values to NoData\n",
    "    data[data == 0] = nodata_value\n",
    "\n",
    "    # Write the modified array back to the DEM\n",
    "    band.WriteArray(data)\n",
    "\n",
    "    # Close the dataset\n",
    "    dem = None\n",
    "\n",
    "def subtract_dems(base_dem_path, aligned_dem_path, subtracted_dem_path):\n",
    "    # Open the base (reference) DEM\n",
    "    base_dem = gdal.Open(base_dem_path)\n",
    "    base_band = base_dem.GetRasterBand(1)\n",
    "    base_array = base_band.ReadAsArray()\n",
    "\n",
    "    # Open the aligned DEM\n",
    "    aligned_dem = gdal.Open(aligned_dem_path)\n",
    "    aligned_band = aligned_dem.GetRasterBand(1)\n",
    "    aligned_array = aligned_band.ReadAsArray()\n",
    "\n",
    "    # Perform subtraction (aligned DEM minus base DEM)\n",
    "    nodata_value = base_band.GetNoDataValue()\n",
    "    difference = np.where((base_array == nodata_value) | (aligned_array == nodata_value), nodata_value, aligned_array - base_array)\n",
    "\n",
    "    # Create a new raster file for the subtracted DEM\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    subtracted_dem = driver.Create(subtracted_dem_path, base_dem.RasterXSize, base_dem.RasterYSize, 1, gdal.GDT_Float32)\n",
    "    subtracted_dem.SetProjection(base_dem.GetProjection())\n",
    "    subtracted_dem.SetGeoTransform(base_dem.GetGeoTransform())\n",
    "\n",
    "    # Write the difference array to the new DEM file\n",
    "    subtracted_band = subtracted_dem.GetRasterBand(1)\n",
    "    subtracted_band.SetNoDataValue(nodata_value)\n",
    "    subtracted_band.WriteArray(difference)\n",
    "\n",
    "    # Close the datasets\n",
    "    base_dem = None\n",
    "    aligned_dem = None\n",
    "    subtracted_dem = None\n",
    "\n",
    "# Specify your DEM file paths\n",
    "base_dem_path = r\"D:\\Bylot\\summer_2017\\Blocs\\2017_ravin_DEM.tif\"\n",
    "target_dem_path = r\"D:\\Bylot\\Summer_2022\\Blocs\\ravin_corrige.tif\"\n",
    "aligned_dem_path = r\"D:\\Bylot\\Summer_2022\\Blocs\\ravin_corrige_aligned4.tif\"\n",
    "subtracted_dem_path = r\"D:\\Bylot\\ravin_difference4.tif\"\n",
    "\n",
    "# Align the DEMs\n",
    "align_dems(base_dem_path, target_dem_path, aligned_dem_path)\n",
    "\n",
    "# Subtract the aligned DEM from the base DEM\n",
    "subtract_dems(base_dem_path, aligned_dem_path, subtracted_dem_path)\n",
    "\n",
    "print(\"DEM alignment, modification, and subtraction complete.\")\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Print computed statistics\n",
    "print(f\"Mean Difference: {mean_difference}\")\n",
    "print(f\"Median Difference: {median_difference}\")\n",
    "print(f\"Standard Deviation: {std_deviation}\")\n",
    "print(f\"Minimum Difference: {min_difference}\")\n",
    "print(f\"Maximum Difference: {max_difference}\")"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Specify the path to the subtracted DEM file\n",
    "subtracted_dem_path = r\"D:\\Bylot\\ravin_difference4.tif\"\n",
    "\n",
    "# Open the DEM with Rasterio\n",
    "with rio.open(subtracted_dem_path) as src:\n",
    "    elevation_difference = src.read(1)\n",
    "    # Set masked or NoData values to np.nan\n",
    "    nodata_value = src.nodata\n",
    "    elevation_difference[elevation_difference == nodata_value] = np.nan\n",
    "\n",
    "# Flatten the array and remove NaN values for statistics\n",
    "flattened_elevation = elevation_difference.flatten()\n",
    "flattened_elevation = flattened_elevation[~np.isnan(flattened_elevation)]\n",
    "\n",
    "# Compute basic statistics\n",
    "mean_difference = np.mean(flattened_elevation)\n",
    "median_difference = np.median(flattened_elevation)\n",
    "std_deviation = np.std(flattened_elevation)\n",
    "min_difference = np.min(flattened_elevation)\n",
    "max_difference = np.max(flattened_elevation)\n",
    "\n",
    "# Round the statistics to two decimal places\n",
    "mean_difference = round(mean_difference, 2)\n",
    "median_difference = round(median_difference, 2)\n",
    "std_deviation = round(std_deviation, 2)\n",
    "min_difference = round(min_difference, 2)\n",
    "max_difference = round(max_difference, 2)\n",
    "\n",
    "# Define the bin range from -0.6 to 0.6 with a bin size of 0.02\n",
    "bin_range = np.arange(-0.6, 0.62, 0.02)\n",
    "\n",
    "# Separate positive and negative elevation differences\n",
    "positive_elevation = flattened_elevation[flattened_elevation > 0]\n",
    "negative_elevation = flattened_elevation[flattened_elevation < 0]\n",
    "\n",
    "# Create histograms for positive and negative values\n",
    "plt.figure(figsize=(12, 8), dpi=300)\n",
    "plt.hist([positive_elevation, negative_elevation], bins=bin_range, color=['blue', 'red'], edgecolor='black', label=['Positive', 'Negative'])\n",
    "plt.title('Histogram of Elevation Differences (-0.6 to 0.6 with 0.02 bin size)')\n",
    "plt.xlabel('Elevation Difference')\n",
    "plt.ylabel('Frequency')\n",
    "plt.xticks(np.arange(-0.6, 0.62, 0.1))\n",
    "plt.legend()\n",
    "plt.show()\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "subtracted_dem_path = r\"D:\\Bylot\\ravin_difference4.tif\"\n",
    "\n",
    "with rasterio.open(subtracted_dem_path) as dem:\n",
    "    dem_data = dem.read(1)  # Read the first band\n",
    "    nodata = dem.nodata  # Get the NoData value\n",
    "\n",
    "# Mask the NoData values\n",
    "dem_data_masked = np.ma.masked_where(dem_data == nodata, dem_data)\n",
    "\n",
    "# Plotting\n",
    "plt.figure(figsize=(20, 14))\n",
    "plt.imshow(dem_data_masked, cmap=\"RdYlGn_r\", vmin=-0.6, vmax=0.6)  # Using original \"coolwarm\" colormap\n",
    "cbar = plt.colorbar(label=\"Delta z (m)\")\n",
    "cbar.ax.tick_params(labelsize=10)\n",
    "cbar.set_label(\"Elevation change (m)\", size=14)# Change font size of color bar ticks\n",
    "plt.title(\"Delta Z from 2022-2017 DEM\", fontsize=22)  # Change title font size\n",
    "#plt.xlabel(fontsize=12)\n",
    "#plt.ylabel(fontsize=12)\n",
    "plt.xticks(fontsize=10)\n",
    "plt.yticks(fontsize=10)\n",
    "plt.savefig(r\"D:\\Bylot\\ravin_dem_difference.jpeg\", dpi=300)"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}