Crop_Landsat7_8_NYC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a2473ebe",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NYC bbox : [563086.4981809029, 4483346.578300821, 609468.3374729961, 4530140.279354534]\n",
      "74\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n",
      "(563086.4981809029, 30.0, 0.0, 4530140.279354534, 0.0, -30.0)\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "Dandan Wei\n",
    "Jan 19, 2022\n",
    "\n",
    "This script crop the Landsat EVI and QA_pixel into the NYC shape, using the Landsat projection.\n",
    "Input dataset:\n",
    "        Landsat7/8 SR_EVI tiff\n",
    "        Landsat7/8 QA_Pixel tiff\n",
    "Output dataset:\n",
    "        Cropped Landsat7/8 SR_EVI tiff\n",
    "        Cropped Landsat7/8 QA_Pixel tiff\n",
    "\n",
    "'''\n",
    "\n",
    "import gdal, ogr, osr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import rasterio\n",
    "import math\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "import glob\n",
    "\n",
    "################################################################\n",
    "#                                                              #\n",
    "#               Define directories and filenames               #\n",
    "#                                                              #\n",
    "################################################################\n",
    "DataDir     = \"/data0/dwei/Landsat/2021/extracted/\" \n",
    "CropDir     = \"/data0/dwei/Landsat/2021/cropped_to_NYC/\" \n",
    "FigDir_temp = \"/data0/dwei/temp/\"\n",
    "ShpeDir     = \"/data0/dwei/NYCshapefiles/\"\n",
    "NYCshapefile= 'geo_export_a8b5d138-820a-4931-8d11-fe8604ea5bd1.shp'\n",
    "\n",
    "################################################################\n",
    "#                                                              #\n",
    "#               Crop the Landsat data into NYC                 #\n",
    "#                                                              #\n",
    "################################################################\n",
    "def get_NYC_bbox(shapefile, target_EPSG):\n",
    "    ''' \n",
    "    To get the bounding box of NYC\n",
    "    '''   \n",
    "    # define shapefile driver\n",
    "    drv         = ogr.GetDriverByName('ESRI Shapefile') \n",
    "    ds_in       = drv.Open(shapefile,0) \n",
    "    lyr_in      = ds_in.GetLayer()            # grab layer\n",
    "    sourceprj   = lyr_in.GetSpatialRef()      # source projection\n",
    "    targetprj   = osr.SpatialReference()\n",
    "    targetprj.ImportFromEPSG(target_EPSG)     # target projection\n",
    "    transform   = osr.CoordinateTransformation(sourceprj, targetprj)\n",
    "\n",
    "    shp         = lyr_in.GetExtent()\n",
    "    ll_x,ll_y,_ = transform.TransformPoint(shp[0],shp[2]) # reprojected bounds for city\n",
    "    ur_x,ur_y,_ = transform.TransformPoint(shp[1],shp[3]) # reprojected bounds for city\n",
    "    bbox        = [ll_x, ll_y, ur_x, ur_y]    # bounding box\n",
    "    # print(\"The original prj is :\\n\", sourceprj)\n",
    "    # print(\"The target   prj is :\\n\", targetprj)\n",
    "    print('NYC bbox :', bbox)\n",
    "    return bbox  \n",
    "\n",
    "def get_geotransform(raster_file_path):\n",
    "    ''' \n",
    "    To print the GeoTransform of a raster file\n",
    "    '''\n",
    "    ds = gdal.Open(raster_file_path)\n",
    "    gt = ds.GetGeoTransform()\n",
    "    print(gt)\n",
    "    return gt \n",
    "\n",
    "# Parameters  \n",
    "target_xRes = 30.0\n",
    "target_yRes = 30.0\n",
    "target_EPSG = 32618    # Landsat projection\n",
    "bbox = get_NYC_bbox(ShpeDir+NYCshapefile, target_EPSG) # NYC minx miny maxx maxy\n",
    "\n",
    "EVI_files = glob.glob(DataDir+'*SR_EVI*')\n",
    "QA_files = glob.glob(DataDir+'*QA_PIXEL*') \n",
    "Band_files = glob.glob(DataDir+'*SR_B*') \n",
    "# print(np.size(QA_files))\n",
    "# print(np.size(QA_files))\n",
    "print(np.size(Band_files))\n",
    "\n",
    "# Crop the Landsat EVI raster files\n",
    "for file in EVI_files:\n",
    "    filename = file[len(DataDir):]\n",
    "    output_evi = CropDir + 'NYC_' + filename\n",
    "#     print(output_evi)\n",
    "    gdal.Warp(output_evi, \\\n",
    "          file, \\\n",
    "          dstSRS = \"EPSG:32618\", \\\n",
    "          xRes = target_xRes, yRes = target_yRes, \\\n",
    "          outputBounds=bbox)  # outputBounds (minX, minY, maxX, maxY)\n",
    "    get_geotransform(output_evi)  \n",
    "    \n",
    "# Crop the Landsat QA_pixel raster files\n",
    "for file in QA_files:\n",
    "    filename = file[len(DataDir):]\n",
    "    output   = CropDir + 'NYC_' + filename\n",
    "#     print(output)\n",
    "    gdal.Warp(output, \\\n",
    "          file, \\\n",
    "          dstSRS = \"EPSG:32618\", \\\n",
    "          xRes = target_xRes, yRes = target_yRes, \\\n",
    "          outputBounds=bbox)  # outputBounds (minX, minY, maxX, maxY)\n",
    "    get_geotransform(output)\n",
    "\n",
    "# Crop the Landsat 7/8 band files\n",
    "for file in Band_files:\n",
    "    filename = file[len(DataDir):]\n",
    "    output   = CropDir + 'NYC_' + filename\n",
    "#     print(output)\n",
    "    gdal.Warp(output, \\\n",
    "          file, \\\n",
    "          dstSRS = \"EPSG:32618\", \\\n",
    "          xRes = target_xRes, yRes = target_yRes, \\\n",
    "          outputBounds=bbox)  # outputBounds (minX, minY, maxX, maxY)\n",
    "    get_geotransform(output)\n",
    "    \n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}