SamplingNoLiqPts_NoLiqEvents.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------------------------------------------------------------------------\n",
    "# fishnet_nonLiq.py\n",
    "# Created on: 08/18/2020\n",
    "# Author: Michelle Meyer\n",
    "# Description: Creates non-liquefaction points for earthquakes with liquefaction\n",
    "#    observances using the fishnet tool. Points are equally spaced 1.5 km apart\n",
    "#    within a 15 km buffer of the liquefaction points but outside of a 1 km buffer.\n",
    "#    NOTE: To sample non-liquefaction points for earthquakes with liquefaction \n",
    "#    observances, see fishnet_nonLiq.py\n",
    "#\n",
    "# ---------------------------------------------------------------------------\n",
    "\n",
    "# Import system modules\n",
    "import arcpy\n",
    "from arcpy import env\n",
    "import pandas as pd\n",
    "import os\n",
    "from arcpy.sa import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change working directory\n",
    "os.chdir(\"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2020_non_liq\\\\z_epicenters\")\n",
    "\n",
    "# Load csv file with folder names\n",
    "eq_list = pd.read_csv(\"non_liq_epicenters.csv\")\n",
    "num_eq = eq_list.shape[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "43.08306566600173\n",
      "13.676319982932643\n",
      "42.36295622760764\n",
      "12.699680017067358\n"
     ]
    }
   ],
   "source": [
    "# How to find extent of each buffer for a fishnet\n",
    "shapeName = arcpy.Describe(forty_km_buffer).shapeFieldName\n",
    "rows = arcpy.SearchCursor(forty_km_buffer)\n",
    "for row in rows:\n",
    "    feat = row.getValue(shapeName)\n",
    "    extent = feat.extent\n",
    "    print(extent.YMax)\n",
    "    print(extent.XMax)\n",
    "    print(extent.YMin)\n",
    "    print(extent.XMin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully completed CentralItaly\n"
     ]
    }
   ],
   "source": [
    "for x in range(0,1):\n",
    "    eq_name = eq_list['NewFolder'][x]\n",
    "\n",
    "    # create path for final shapefile and intermediate files \n",
    "    output_path = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2021_liq_Alex\\\\\" + eq_name + \"\\\\\"\n",
    "    #os.mkdir(output_path)    \n",
    "    #temp_path = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2021_liq_Alex\\\\\" + eq_name + \"\\\\\"\n",
    "    #os.mkdir(temp_path)\n",
    "    \n",
    "    # Grab epicenter  \n",
    "    epicenter = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2021_liq_Alex\\\\z_epicenters\\\\\" + eq_name + \".shp\"\n",
    "\n",
    "    # find spatial reference\n",
    "    raster_reference = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\global_data_2020\\\\GMTED\\\\\" + eq_name + \"\\\\\" + eq_name + \"_DEM_proj.tif\"\n",
    "    spatial_ref = arcpy.Describe(raster_reference).spatialReference\n",
    "    \n",
    "    out_epicenter = output_path + eq_name + \"_epi_proj.shp\"\n",
    "    arcpy.management.Project(epicenter, out_epicenter, spatial_ref)\n",
    "    \n",
    "    forty_km_buffer = temp_path + \"\\\\40km_buffer.shp\"\n",
    "    arcpy.Buffer_analysis(epicenter, forty_km_buffer, \"40000 meters\", \"\", \"\", \"ALL\")\n",
    "    \n",
    "    print(\"Successfully completed \"+eq_name)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in range(0,num_eq):\n",
    "    eq_name = eq_list['NewFolder'][x]\n",
    "\n",
    "    # create path for final shapefile and intermediate files \n",
    "    output_path = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2021_liq_Alex\\\\\" + eq_name + \"\\\\\"\n",
    "    os.mkdir(output_path)    \n",
    "    temp_path = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2021_liq_Alex\\\\\" + eq_name + \"\\\\\"\n",
    "    os.mkdir(temp_path)\n",
    "    \n",
    "    # Grab epicenter  \n",
    "    epicenter = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\Liquefaction_inventories\\\\2020_non_liq\\\\z_epicenters\\\\\" + eq_name + \".shp\"\n",
    "\n",
    "    # find spatial reference\n",
    "    raster_reference = \"R:\\\\Liquefaction_RemoteSensing\\\\Projects\\\\GeospatialModel\\\\global_data_2020\\\\GMTED\\\\\" + eq_name + \"\\\\\" + eq_name + \"_DEM_proj.tif\"\n",
    "    spatial_ref = arcpy.Describe(raster_reference).spatialReference\n",
    "    \n",
    "    \n",
    "    arcpy.management.Project(epicenter, out_dataset, spatial_ref)\n",
    "    forty_km_buffer = temp_path + \"\\\\40km_buffer.shp\"\n",
    "    arcpy.Buffer_analysis(epicenter, forty_km_buffer, \"40000 meters\", \"\", \"\", \"ALL\")\n",
    "    \n",
    "    \n",
    "    env.workspace = temp_path\n",
    "    env.outputCoordinateSystem = spatial_ref\n",
    "    \n",
    "    \n",
    "    \n",
    "    left_corner = arcpy.GetRasterProperties_management(forty_km_buffer, \"LEFT\")\n",
    "    right_corner = arcpy.GetRasterProperties_management(forty_km_buffer, \"RIGHT\")\n",
    "    bottom_corner = arcpy.GetRasterProperties_management(forty_km_buffer, \"BOTTOM\")\n",
    "    top_corner = arcpy.GetRasterProperties_management(forty_km_buffer, \"TOP\")\n",
    "      \n",
    "    \n",
    "    # create fishnet\n",
    "    originCoordinate = left_corner[0] + ' ' + bottom_corner[0]\n",
    "    yAxisCoordinate = left_corner[0] + ' ' + str(float(bottom_corner[0])+10)\n",
    "     \n",
    "    cellSizeWidth = '1500'\n",
    "    cellSizeHeight = '1500'\n",
    "    \n",
    "    numRows =  '0'\n",
    "    numColumns = '0'\n",
    "    \n",
    "    oppositeCorner = right_corner[0] + ' ' + top_corner[0]\n",
    "    \n",
    "    fishnet_out = temp_path + \"\\\\fishnet.shp\"\n",
    "    fishnet_points = temp_path + \"\\\\fishnet_label.shp\"\n",
    "    arcpy.CreateFishnet_management(fishnet_out, originCoordinate, yAxisCoordinate, cellSizeWidth, cellSizeHeight, numRows, numColumns, oppositeCorner, 'LABELS', '#', 'POLYGON')\n",
    "   \n",
    "    # clip to sample area shape (created from buffer tool)\n",
    "    fishnet_clipped = temp_path + \"\\\\fishnet_clipped.shp\"\n",
    "    arcpy.Clip_analysis(fishnet_points, forty_km_buffer, fishnet_clipped)\n",
    "    \n",
    "    # change points to GCS\n",
    "    out_coordinate_system = arcpy.SpatialReference('WGS 1984')\n",
    "    non_liq_points = output_path + \"\\\\non_liq_point.shp\"\n",
    "    arcpy.Project_management(fishnet_clipped, non_liq_points, out_coordinate_system)\n",
    "    \n",
    "    # add X and Y coordinates\n",
    "    env.outputCoordinateSystem = out_coordinate_system\n",
    "    arcpy.AddXY_management(non_liq_points)\n"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}