generate_CEHZ_variables.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import modules\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import rasterstats\n",
    "import rasterio\n",
    "from shapely.geometry import Polygon, LineString, Point \n",
    "import numpy as np \n",
    "import tqdm\n",
    "from tqdm.notebook import tqdm_notebook\n",
    "import time\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "warnings.simplefilter(action='ignore', category=UserWarning) \n",
    "pd.options.mode.chained_assignment = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### STEP 1 - Find nearest wave point for each transect ####\n",
    "# read transect points to gpd and tidy up\n",
    "#set engine to pyogrio\n",
    "engine='pyogrio'\n",
    "transects = gpd.read_file('data/Intersect_v2.shp', engine=engine).to_crs(2193)\n",
    "\n",
    "# list of attributes to keep\n",
    "to_keep = ['BaselineID', 'TransectID', 'NSM', 'EPR', 'WLR', 'WCI90', 'DSASdate', 'geometry']\n",
    "# iterate over cols and drop unwanted cols\n",
    "to_drop = []\n",
    "for col in transects.columns:\n",
    "    if col not in to_keep:\n",
    "        to_drop.append(col)\n",
    "print(to_drop)\n",
    "transects.drop(to_drop, axis=1, inplace=True)\n",
    "\n",
    "print(transects.value_counts('TransectID'))\n",
    "transects.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read wave data incl CD as gpd dataframe\n",
    "wave_points_df = pd.read_csv('data/wave_points_nztm.csv')\n",
    "wave_points_gdf = gpd.GeoDataFrame(wave_points_df, geometry=gpd.points_from_xy(wave_points_df.X, wave_points_df.Y), crs=2193)\n",
    "wave_points_gdf.rename(columns={'X': 'wave_X', 'Y': 'wave_Y', 'ID': 'waveID'}, inplace=True)\n",
    "wave_points_gdf.head()\n",
    "\n",
    "# Perform nearest sjoin to find nearest wave_point for each transect\n",
    "transect_wp = gpd.sjoin_nearest(transects, wave_points_gdf, how='inner')\n",
    "transect_wp.head()\n",
    "\n",
    "### END OF STEP 1 ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### STEP 2 - Create .shp for transects from intersection to wave point ####\n",
    "# create new gdf where geometry is Linestring from transect intersect to wave point\n",
    "# function to generate Linestring geometry\n",
    "def get_transect_geom(row):\n",
    "    eov_x = row.geometry.x\n",
    "    eov_y = row.geometry.y\n",
    "    wave_x = row.wave_X\n",
    "    wave_y = row.wave_Y\n",
    "    return  LineString([(eov_x, eov_y), (wave_x, wave_y)])\n",
    "\n",
    "def get_eov_x(row):\n",
    "    eov_x = row.geometry.x\n",
    "    return  eov_x\n",
    "\n",
    "def get_eov_y(row):\n",
    "    eov_y = row.geometry.y\n",
    "    return  eov_y\n",
    "\n",
    "# create new column containing Linestring from transect to wave point\n",
    "transect_wp['transect_wp_geom'] = transect_wp.apply(lambda row : get_transect_geom(row), axis=1)\n",
    "print(transect_wp.head())\n",
    "\n",
    "# store intersect \n",
    "transect_wp['intsect_X'] = transect_wp.apply(lambda row : get_eov_x(row), axis=1)\n",
    "transect_wp['intsect_Y'] = transect_wp.apply(lambda row : get_eov_y(row), axis=1)\n",
    "\n",
    "# set geometry to transect_wp Linestring\n",
    "# tidy up gdf and set geometry to transect\n",
    "transect_wp.set_geometry('transect_wp_geom', inplace=True, crs=2193)\n",
    "#transect_wp.drop(['index_right'], axis=1, inplace=True)\n",
    "\n",
    "transect_wp.drop(columns=['geometry'], axis=1, inplace=True)\n",
    "transect_wp.head()\n",
    "\n",
    "### END OF STEP 2 ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### STEP 3 FUNCTIONS ### \n",
    "# function to find depth corresponding to closure depth\n",
    "def get_depth(df):\n",
    "    # define cd\n",
    "    cd = df.CD\n",
    "    \n",
    "    # return the bathy_dpth value closest to CD \n",
    "    bthy_dpth = df.iloc[(df['bathy_dpth']--abs(cd.iloc[0])).abs().argsort()[:1]] # using -abs() to ensure cd depth values match bathy_dpth\n",
    "    \n",
    "    return bthy_dpth\n",
    "\n",
    "# function to find length between\n",
    "def get_length(df):\n",
    "    # create gdf to find length of transect from elevation to depth\n",
    "    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.intsect_X, df.intsect_Y), crs=2193)\n",
    "   \n",
    "    # redefine geometry from bathy_depth to eov_point and set as geometry\n",
    "    gdf['transect'] = LineString([(df.intsect_X, df.intsect_Y), (df.coord_X, df.coord_Y)])\n",
    "    gdf.set_geometry('transect', inplace=True, crs=2193)\n",
    "    \n",
    "    # calc length of transect\n",
    "    gdf['length'] = gdf.length\n",
    "    \n",
    "    return gdf\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### STEP 3 - FIND DEPTH VALUES ####\n",
    "\n",
    "### STEP 3 FUNCTIONS ### \n",
    "# function to find depth corresponding to closure depth\n",
    "def get_depth(df):\n",
    "    # define cd\n",
    "    cd = df.CD\n",
    "    \n",
    "    # return the bathy_dpth value closest to CD \n",
    "    bthy_dpth = df.iloc[(df['bathy_dpth']--abs(cd.iloc[0])).abs().argsort()[:1]] # using -abs() to ensure cd depth values match bathy_dpth\n",
    "    \n",
    "    return bthy_dpth\n",
    "\n",
    "# function to find length between\n",
    "def get_length(df):\n",
    "    # create gdf to find length of transect from elevation to depth\n",
    "    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.intsect_X, df.intsect_Y), crs=2193)\n",
    "   \n",
    "    # redefine geometry from bathy_depth to eov_point and set as geometry\n",
    "    gdf['transect'] = LineString([(df.intsect_X, df.intsect_Y), (df.coord_X, df.coord_Y)])\n",
    "    gdf.set_geometry('transect', inplace=True, crs=2193)\n",
    "    \n",
    "    # calc length of transect\n",
    "    gdf['length'] = gdf.length\n",
    "    \n",
    "    return gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIND BATHY VALUES WITHIN 250m OF EACH TRANSECT\n",
    "# read bathymetry centroids to gdf\n",
    "start = time.time()\n",
    "bathy = gpd.read_file('data/Bathy250m_centroids_nztm.shp', engine=engine).to_crs(2193)\n",
    "end = time.time()\n",
    "print('time to read bathmetry: ', end - start, 'seconds')\n",
    "\n",
    "start = time.time()\n",
    "# get coordinates for bathy depth for next step\n",
    "bathy['coord_X'] = bathy.apply(lambda row : get_eov_x(row), axis=1)\n",
    "bathy['coord_Y'] = bathy.apply(lambda row : get_eov_y(row), axis=1)\n",
    "\n",
    "# generate buffers for each transect and set_geometry to buffer polygons\n",
    "transect_wp['buffer'] = transect_wp.buffer(250)\n",
    "transect_wp.set_geometry('buffer', inplace=True, crs=2193)\n",
    "# transect_wp.drop(['index_right'], axis=1, inplace=True)\n",
    "transect_wp.head()\n",
    "\n",
    "# spatial join to find all bathy_points that lie within transect buffers \n",
    "points_within = gpd.sjoin(bathy, transect_wp,  how='inner', predicate='within')\n",
    "# drop join columns\n",
    "points_within.drop(['index_right'], axis=1, inplace=True)\n",
    "\n",
    "end = time.time()\n",
    "print('time to perform spatial join: ', end - start, 'seconds')\n",
    "points_within.head()\n",
    "\n",
    "\n",
    "### RETURN BATHY DEPTH and length of each transect \n",
    "# iterate over each transect and find the bathy_dpth closest to cd and eov elevation\n",
    "# create empty list of df with same columns\n",
    "df_list = []\n",
    "# get list of transect IDs\n",
    "id_list = points_within.TransectID.unique().tolist()\n",
    "\n",
    "# gen start time\n",
    "start = time.time()\n",
    "\n",
    "for i in tqdm_notebook(id_list, desc='Finding bathymetry depth: '):\n",
    "    # return rows corresponding to a given transect\n",
    "    rows = points_within.query('TransectID == @i')\n",
    "\n",
    "    # get bathymetry depth\n",
    "    depth = get_depth(rows)\n",
    "\n",
    "    # get length \n",
    "    gdf = get_length(depth)\n",
    "\n",
    "    # append to dataframe list\n",
    "    df_list.append(gdf)\n",
    "\n",
    "end = time.time()\n",
    "elapsed = end - start\n",
    "print('run time: ', elapsed, 'seconds')\n",
    "\n",
    "# concat dataframe list and write to csv file\n",
    "df_bthy_dpth = pd.DataFrame(pd.concat(df_list))  \n",
    "\n",
    "# save to csv\n",
    "df_bthy_dpth.to_csv('all_depth_length.csv')\n",
    "df_bthy_dpth.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### STEP 4 - CALC ELEVATION BASED ON EOV POINT\n",
    "\n",
    "## FUNCTIONS ##\n",
    "# function to find EOV elevation from LIDAR\n",
    "def get_elevation(df_row, raster):\n",
    "    # extract xy from columns\n",
    "    x = df_row.intsect_X\n",
    "    y = df_row.intsect_Y\n",
    "    # create point from x, y\n",
    "    point = Point(x,y)\n",
    "    # point query from rasterstats returning elevation value\n",
    "    ele_val = rasterstats.point_query(point, raster)\n",
    "    print(df_row.TransectID)\n",
    "    print(\"Raster value on point %.2f \\n\"%ele_val[0])\n",
    "    \n",
    "    return ele_val[0]\n",
    "\n",
    "## END OF FUNCTIONS ##\n",
    "\n",
    "# read in and join LIDAR tile ids\n",
    "tile_id = gpd.read_file('/path/to/lidar/tiles/shapefile').to_crs(2193)\n",
    "\n",
    "# spatial join to get tile_id\n",
    "# spatial join to find all bathy_points that lie within transect buffers \n",
    "elevation_points = gpd.sjoin(df_bthy_dpth, tile_id,  how='inner', predicate='within')\n",
    "\n",
    "elevation_points.drop(['field_1','index_right', 'topo50_cov'], axis=1, inplace=True)\n",
    "\n",
    "print(elevation_points.value_counts('sheet_code'))\n",
    "\n",
    "elevation_points.head()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create list of tile_id to iterate over\n",
    "tile_id_list = elevation_points.sheet_code.unique().tolist()\n",
    "\n",
    "# create blank list to save df to\n",
    "df_list = []\n",
    "\n",
    "start = time.time()\n",
    "for i in  tqdm_notebook(tile_id_list,desc='Finding elevation: '):  \n",
    "    rows = elevation_points.query('sheet_code == @i')\n",
    "    # read LUT containing LiDAR tile file paths\n",
    "    tile_lut = pd.read_csv('data/LIDARTILES.csv')\n",
    "    # find tile using LUT of tile file paths\n",
    "    tile_id = rows.sheet_code\n",
    "    #print('finding tile: ' + tile_id)\n",
    "    lidar_fp = ''\n",
    "    for i in range(len(tile_lut)):\n",
    "           if tile_lut.loc[i, \"TILES\"] == tile_id.iloc[0]:\n",
    "                 lidar_fp += tile_lut.loc[i, \"FP\"]                \n",
    "    print('opening LiDAR tile: ' + lidar_fp)\n",
    "    rows['elevation'] = rows.apply(lambda row: get_elevation(row, lidar_fp), axis=1)\n",
    "    df_list.append(rows)\n",
    "\n",
    "end = time.time()\n",
    "print('time to find elevation values from lidar: ', end - start, 'seconds') \n",
    "\n",
    "\n",
    "# # concat dataframe list and write to csv df_list))\n",
    "df_nearly_there = pd.DataFrame(pd.concat(df_list))\n",
    "df_nearly_there.to_csv('all_CEHZ_variables.csv')\n",
    "df_nearly_there.head()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}