rasterSampler.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e7602d3c-218f-4315-9d71-587a314676f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import rasterio\n",
    "import shapely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "76abcfe8-6cb8-4382-a391-9b7aa1609b5d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████████████████████████| 2456/2456 [00:02<00:00, 823.85it/s]\n",
      "100%|██████████████████████████████████████████████████████████████████████████████████████████| 6022112/6022112 [00:19<00:00, 305682.63it/s]\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(\"field.tif\") as src:\n",
    "    data = src.read()\n",
    "    bounds = src.bounds\n",
    "    crs = src.crs\n",
    "\n",
    "bands, height_in_pixels, width_in_pixels = data.shape\n",
    "\n",
    "pixel_height = np.abs(bounds.top-bounds.bottom)/height_in_pixels\n",
    "pixel_width = np.abs(bounds.left-bounds.right)/width_in_pixels\n",
    "\n",
    "pixel_centroids = []\n",
    "for h in tqdm(range(height_in_pixels),position=0,leave=True):\n",
    "    for w in range(width_in_pixels):\n",
    "        pixel_centroids.append(\n",
    "            [bounds.left+(w+1/2)*pixel_width,\n",
    "             bounds.bottom+(h+1/2)*pixel_height]\n",
    "        )\n",
    "        \n",
    "pixel_centroids = gpd.GeoSeries([shapely.geometry.Point(e) for e in tqdm(pixel_centroids)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c89e9d0d-fb87-45b0-b5dd-418f5c6b6dae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get pixel centroids within cutline\n",
    "masks = []\n",
    "for geom in gpd.read_file(\"cutline.geojson\").geometry:\n",
    "    masks.append(\n",
    "        pixel_centroids.within(geom)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7dc36d27-4e6a-44fd-9814-1c5e854a012a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample the raster for those centroids\n",
    "dfsamples = []\n",
    "for mask in masks:\n",
    "    samples = gpd.GeoDataFrame(geometry = pixel_centroids[mask]).reset_index(drop=True)\n",
    "    with rasterio.open(\"field.tif\") as src:\n",
    "        samples = samples.assign(pixelvals = samples.geometry.apply(lambda p: next(src.sample(list(zip(*p.xy))))))\n",
    "    dfsamples.append(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5cb2c8ce-3c34-40e1-b7a3-01a2659fcf7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# join geometry and sample columns, write coordinate to two columns\n",
    "for i, df in enumerate(dfsamples):\n",
    "    \n",
    "    banddf = pd.DataFrame(df.pixelvals.tolist())\n",
    "    banddf.columns = [f\"band{e+1:02}\" for e in banddf.columns]\n",
    "    \n",
    "    df = df.join(\n",
    "        pd.DataFrame(\n",
    "            df.geometry.map(lambda p: list(*zip(*p.xy))).tolist(),\n",
    "            columns=['x','y']\n",
    "        )\n",
    "    ).join(banddf)\n",
    "    \n",
    "    df = df[[c for c in df.columns if c not in [\"pixelvals\",'geometry']]]\n",
    "    df.to_csv(f\"samples{i}.csv\",index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev311_2",
   "language": "python",
   "name": "dev311_2"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}