lai_pro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterio.crs import CRS\n",
    "from rasterio.enums import Resampling\n",
    "from rasterio.warp import calculate_default_transform, reproject\n",
    "from rasterio.mask import mask\n",
    "from rasterio.features import shapes\n",
    "from shapely.geometry import shape, mapping\n",
    "from shapely.ops import unary_union\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_nc_to_tif(nc_input_path, temp_tif_output_path, ref_raster_path, aligned_tif_path, final_output_path):\n",
    "    # Open the NetCDF file\n",
    "    with rasterio.open(nc_input_path) as nc_file:\n",
    "        meta = nc_file.meta.copy()\n",
    "        meta.update({\n",
    "            'driver': 'GTiff',\n",
    "            'height': nc_file.height,\n",
    "            'width': nc_file.width,\n",
    "            'crs': CRS.from_epsg(4326),\n",
    "            'transform': rasterio.transform.from_origin(nc_file.bounds.left, nc_file.bounds.top, nc_file.res[0], nc_file.res[1])\n",
    "        })\n",
    "\n",
    "        # Read data\n",
    "        data = nc_file.read(1)\n",
    "\n",
    "        # Write data to GeoTIFF\n",
    "        with rasterio.open(temp_tif_output_path, 'w', **meta) as dst:\n",
    "            dst.write(data, 1)\n",
    "\n",
    "    # Align and reproject the tif\n",
    "    with rasterio.open(temp_tif_output_path) as src, rasterio.open(ref_raster_path) as ref:\n",
    "        transform, width, height = calculate_default_transform(\n",
    "            src.crs, ref.crs, ref.width, ref.height, *ref.bounds)\n",
    "        meta = src.meta.copy()\n",
    "        meta.update({\n",
    "            'crs': ref.crs,\n",
    "            'transform': transform,\n",
    "            'width': width,\n",
    "            'height': height,\n",
    "            'driver': 'GTiff'\n",
    "        })\n",
    "        \n",
    "        with rasterio.open(aligned_tif_path, 'w', **meta) as dst:\n",
    "            for i in range(1, src.count + 1):\n",
    "                reproject(\n",
    "                    source=rasterio.band(src, i),\n",
    "                    destination=rasterio.band(dst, i),\n",
    "                    src_transform=src.transform,\n",
    "                    src_crs=src.crs,\n",
    "                    dst_transform=transform,\n",
    "                    dst_crs=ref.crs,\n",
    "                    resampling=Resampling.nearest)\n",
    "\n",
    "    # Clip the raster using the reference raster as a mask\n",
    "    polygon = create_mask_polygon(ref_raster_path)\n",
    "    if polygon:\n",
    "        clip_raster_with_polygon(aligned_tif_path, final_output_path, polygon)\n",
    "\n",
    "def create_mask_polygon(raster_path):\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        image = src.read(1)\n",
    "        mask = image > 0\n",
    "        transform = src.transform\n",
    "        shape_gen = shapes(image, mask=mask, transform=transform)\n",
    "        polygons = [shape(geom) for geom, value in shape_gen if value]\n",
    "        if polygons:\n",
    "            return unary_union(polygons)\n",
    "    return None\n",
    "\n",
    "def clip_raster_with_polygon(raster_path, output_path, polygon):\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        out_image, out_transform = mask(src, [mapping(polygon)], crop=True)\n",
    "        out_meta = src.meta.copy()\n",
    "        out_meta.update({\n",
    "            \"driver\": \"GTiff\",\n",
    "            \"height\": out_image.shape[1],\n",
    "            \"width\": out_image.shape[2],\n",
    "            \"transform\": out_transform\n",
    "        })\n",
    "\n",
    "        with rasterio.open(output_path, 'w', **out_meta) as out_raster:\n",
    "            out_raster.write(out_image)\n",
    "\n",
    "# Define the data types and process each\n",
    "base_path = 'data/aggregated_over_years_model_input/lai_min.nc'\n",
    "ref_raster_path = 'data/dead/deadwood.tif'\n",
    "\n",
    "nc_input = base_path\n",
    "temp_tif_output = 'data/temp/lai_min_temp.tif'\n",
    "aligned_tif_output = 'data/temp/lai_min_aligned.tif'\n",
    "final_output_path = 'data/aggregated_over_years_model_input/lai_min_pro.tif'\n",
    "process_nc_to_tif(nc_input, temp_tif_output, ref_raster_path, aligned_tif_output, final_output_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_path = 'data/aggregated_over_years_model_input/lai_sum_mean.nc'\n",
    "ref_raster_path = 'data/dead/deadwood.tif'\n",
    "\n",
    "nc_input = base_path\n",
    "temp_tif_output = 'data/temp/lai_sum_mean.tif'\n",
    "aligned_tif_output = 'data/temp/lai_sum_mean_aligned.tif'\n",
    "final_output_path = 'data/aggregated_over_years_model_input/lai_sum_mean_pro.tif'\n",
    "process_nc_to_tif(nc_input, temp_tif_output, ref_raster_path, aligned_tif_output, final_output_path)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "babadormann",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}