{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### Microsoft Azure open data program https://azure.microsoft.com/en-us/services/open-datasets/catalog/naip/" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import geopandas as gpd\n", "import os, shutil\n", "import rasterio\n", "import urllib\n", "import re\n", "import fiona.transform\n", "import requests\n", "#from rasterstats import zonal_stats\n", "import numpy as np\n", "import gdal \n", "gdal.SetConfigOption(\"GDAL_HTTP_UNSAFESSL\", \"YES\")\n", "gdal.VSICurlClearCache()" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# Put the quad and aoi files in the infilepath directory, define a different outpath if desired\n", "infilepath = r'Z:\\tiger\\scripts\\FireLine_ImageProcessing\\NAIP'\n", "outpath = infilepath\n", "\n", "#text to append to each NAIP quad's complete analysis\n", "outfilen_append = '_msavi.tif'\n", "\n", "# shapefile from microsoft Azure with NAIP quads\n", "# https://naipeuwest.blob.core.windows.net/naip/v002/ca/2018/ca_shpfl_2018/index.html\n", "inNAIPquad = 'NAIP_18_CA.shp'\n", "\n", "#file you create with polygons/points of interest. If using the whole state: aoi_file = 'centroid'\n", "aoi_file = 'sample_aoi.gpkg'" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "# NAIP metadata\n", "year = '2018' \n", "state = 'ca' #lowercase\n", "resolution = '060cm'\n", "blob_root = 'https://naipeuwest.blob.core.windows.net/naip'" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "## Read the aoi file. if using the whole state, the code will create a list of tiles based on centroids later\n", "if aoi_file == 'centroid':\n", " aoi = 'centroid'\n", "else:\n", " aoi = gpd.read_file(os.path.join(infilepath, aoi_file))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "def createListofTiles(inpath,NAIPquad,inaoi):\n", " # NAIP Quad file\n", " ntiles = gpd.GeoDataFrame.from_file(os.path.join(inpath,NAIPquad))\n", " #if inaoi == 'centroid':\n", " if isinstance(inaoi, str):\n", " print('centroid')\n", " aoicreate = gpd.GeoDataFrame(ntiles.geometry.centroid, columns=['geometry'])\n", " inaoi = gpd.sjoin(aoicreate, ntiles, how=\"inner\", op='intersects')\n", " inaoi = inaoi[['geometry']] \n", " \n", " if inaoi.crs != ntiles.crs:\n", " print(inaoi.crs,ntiles.crs)\n", " c = ntiles.crs.dtypes\n", " print(c)\n", " inaoi= inaoi.to_crs(c)\n", " \n", " \n", " # sort so that the images with the most aoi's get processed first, if applicable\n", " if str((aoi.geom_type == 'Point')[0]) == 'True':\n", " # Intersection with AOI to the names of all the tiles needed\n", " quad_intersection = gpd.sjoin(ntiles, inaoi, how=\"right\", op='intersects')\n", " print(quad_intersection)\n", " \n", " print('Sorting by counts of points in each quad')\n", " quad_intersection['counts'] = quad_intersection.groupby(['FileName'])['FileName'].transform('count')\n", " quad_intersection = quad_intersection.sort_values(by='counts',ascending=False)\n", " else:\n", " # Intersection with AOI to the names of all the tiles needed\n", " quad_intersection = gpd.sjoin(inaoi, ntiles, how=\"left\", op='intersects')\n", " print(quad_intersection)\n", " \n", " quadlist = quad_intersection.FileName.unique()\n", " return quadlist" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " field geometry index_right \\\n", "0 sample2 POLYGON ((-13533555.220 4763989.707, -13534804... 2964 \n", "0 sample2 POLYGON ((-13533555.220 4763989.707, -13534804... 2963 \n", "\n", " AREA PERIMETER XCOORD YCOORD ST QQNAME QKEY ... \\\n", "0 0.004 0.25 121.5000 39.25 CA HONCUT SE 3915001213000 ... \n", "0 0.004 0.25 121.5625 39.25 CA HONCUT SW 3915001213345 ... \n", "\n", " QKEY_1 BAND USGSID QUAD UTM RES SrcImgDate VerDate \\\n", "0 3.915001e+12 M4B 3912144 SE 10 60 20180718 20190210 \n", "0 3.915001e+12 M4B 3912144 SW 10 60 20180721 20190210 \n", "\n", " FileName CELL_1 \n", "0 m_3912144_se_10_060_20180718_20190210.tif SE3912144 \n", "0 m_3912144_sw_10_060_20180721_20190210.tif SW3912144 \n", "\n", "[2 rows x 33 columns]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "c:\\python27v2\\lib\\site-packages\\numpy\\lib\\function_base.py:2167: RuntimeWarning: invalid value encountered in find_intersects (vectorized)\n", " outputs = ufunc(*inputs)\n" ] } ], "source": [ "ql = createListofTiles(infilepath,inNAIPquad,aoi)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "def calcNDVI(src):\n", " # read red and nir bands, change to dtype to float64 for output ndvi calc\n", " red,nir = src.read(1).astype('float64'),src.read(4).astype('float64')\n", " print('Calculate NDVI')\n", " # Calculate NDVI\n", " ndvi = np.where((nir+red) == 0, 0, (nir-red)/(nir+red))\n", "\n", " # Return reclassified NDVI into fuel (1) or no fuel (0)\n", " ndvi[ndvi>.1] = 1\n", " ndvi[ndvi<=.1] = 0\n", " return ndvi\n", " \n", " \n", "def calcMSAVI2(src):\n", " from math import sqrt\n", " # read red and nir bands, change to dtype to float64 for output calc\n", " red,nir = src.read(1).astype('float64'),src.read(4).astype('float64')\n", " print('Calculate MSAVI2')\n", " msavi2 = (2 * (nir.astype(float) + 1) - np.sqrt((2 * nir.astype(float) + 1)**2 - 8 * (nir.astype(float) - red.astype(float))))/2\n", "\n", " return msavi2\n", " \n", " " ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 of 2\n", "https://naipeuwest.blob.core.windows.net/naip/v002/ca/2018/ca_060cm_2018/39121/m_3912144_se_10_060_20180718_20190210.tif\n", "Calculate MSAVI2\n", "m_3912144_se_10_060_20180718_20190210_msavi.tif Complete\n", "2 of 2\n", "https://naipeuwest.blob.core.windows.net/naip/v002/ca/2018/ca_060cm_2018/39121/m_3912144_sw_10_060_20180721_20190210.tif\n", "Calculate MSAVI2\n", "m_3912144_sw_10_060_20180721_20190210_msavi.tif Complete\n" ] } ], "source": [ "n = 0\n", "\n", "for filen in (ql):\n", " n=n+1\n", " print(str(n) + ' of ' + str(len(ql)))\n", " \n", " # Retrieve the name of the quad file and define link to file\n", " quadrangle = filen[2:7]\n", " filename = filen\n", " tile_url = blob_root + '/v002/' + state + '/' + year + '/' + state + '_' + resolution + \\\n", " '_' + year + '/' + quadrangle + '/' + filename\n", " \n", " print(tile_url)\n", " \n", " with rasterio.open(tile_url) as f:\n", " \n", " #### Update analysis to be run\n", " analysis = calcMSAVI2(f)\n", "\n", " # set transform\n", " affine= f.transform\n", "\n", " outputanalysis = os.path.join(outpath,filename[:-4]+outfilen_append)\n", " with rasterio.open(outputanalysis, 'w', \n", " driver='GTiff',\n", " height=f.shape[0],\n", " width=f.shape[1],\n", " compress='lzw', #lossless\n", " count=1,\n", " dtype='float64', # 'uint8' for scaled ndvi or reclass\n", " crs=f.crs,\n", " transform=affine) as dst:\n", " dst.write(analysis,1)\n", " print(filename[:-4]+ outfilen_append + ' Complete')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 2 }