API tools faq
paste
Guest User

Untitled

a guest
Dec 12th, 2018
50
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 14.16 KB | None | 0 0
raw download clone embed print report
  1. {
  2. "cells": [
  3. {
  4. "cell_type": "markdown",
  5. "metadata": {},
  6. "source": [
  7. "The aim of this recommended practice is to generate a burn severity map for the assessment of the areas affected by the wildfires using the differenced Normalized Burn Ratio (dNBR), which is the difference of the post-fire NBR and the pre-fire NBR. "
  8. ]
  9. },
  10. {
  11. "cell_type": "markdown",
  12. "metadata": {},
  13. "source": [
  14. "Import the required packages"
  15. ]
  16. },
  17. {
  18. "cell_type": "code",
  19. "execution_count": null,
  20. "metadata": {},
  21. "outputs": [],
  22. "source": [
  23. "import osr\n",
  24. "import ogr\n",
  25. "import gdal\n",
  26. "import numpy as np\n",
  27. "import matplotlib\n",
  28. "import matplotlib.pyplot as plt\n",
  29. "import rasterio\n",
  30. "import glob\n",
  31. "from rasterio.plot import show\n",
  32. "from rasterio.mask import mask\n",
  33. "from shapely.geometry import mapping\n",
  34. "import geopandas as gpd\n",
  35. "import math\n"
  36. ]
  37. },
  38. {
  39. "cell_type": "markdown",
  40. "metadata": {},
  41. "source": [
  42. "Function definitions"
  43. ]
  44. },
  45. {
  46. "cell_type": "code",
  47. "execution_count": null,
  48. "metadata": {},
  49. "outputs": [],
  50. "source": [
  51. "def read_band_image(band, path):\n",
  52. " \"\"\"\n",
  53. " This function takes as input the Sentinel-2 band name and the path of the \n",
  54. " folder that the images are stored, reads the image and returns the data as\n",
  55. " an array\n",
  56. " input: band string Sentinel-2 band name\n",
  57. " path string path of the folder\n",
  58. " output: data array (n x m) array of the band image\n",
  59. " spatialRef string projection \n",
  60. " geoTransform tuple affine transformation coefficients\n",
  61. " targetprj spatial reference\n",
  62. " \"\"\"\n",
  63. " a = path+'*B'+band+'*.jp2'\n",
  64. " img = gdal.Open(glob.glob(a)[0])\n",
  65. " data = np.array(img.GetRasterBand(1).ReadAsArray())\n",
  66. " spatialRef = img.GetProjection()\n",
  67. " geoTransform = img.GetGeoTransform()\n",
  68. " targetprj = osr.SpatialReference(wkt = img.GetProjection())\n",
  69. " return data, spatialRef, geoTransform, targetprj\n",
  70. "\n",
  71. "def nbr(band1, band2):\n",
  72. " \"\"\"\n",
  73. " This function takes an input the arrays of the bands from the read_band_image\n",
  74. " function and returns the Normalized Burn ratio (NBR)\n",
  75. " input: band1 array (n x m) array of first band image e.g B8A\n",
  76. " band2 array (n x m) array of second band image e.g. B12\n",
  77. " output: nbr array (n x m) normalized burn ratio\n",
  78. " \"\"\"\n",
  79. " nbr = (band1 - band2) / (band1 + band2)\n",
  80. " return nbr\n",
  81. "\n",
  82. "def dnbr(nbr1,nbr2):\n",
  83. " \"\"\"\n",
  84. " This function takes as input the pre- and post-fire NBR and returns the dNBR\n",
  85. " input: nbr1 array (n x m) pre-fire NBR\n",
  86. " nbr2 array (n x m) post-fire NBR\n",
  87. " output: dnbr array (n x m) dNBR\n",
  88. " \"\"\"\n",
  89. " dnbr = nbr1 - nbr2\n",
  90. " return dnbr\n",
  91. "\n",
  92. "def reproject_shp_gdal(infile, outfile, targetprj):\n",
  93. " \"\"\"\n",
  94. " This function takes as input the input and output file names and the projection\n",
  95. " in which the input file will be reprojected and reprojects the input file using\n",
  96. " gdal\n",
  97. " input: infile string input filename\n",
  98. " outfile string output filename\n",
  99. " targetprj projection (output of function read_band_image)\n",
  100. " \"\"\"\n",
  101. " ## reprojection with gdal \n",
  102. " \n",
  103. " driver = ogr.GetDriverByName(\"ESRI Shapefile\") \n",
  104. " dataSource = driver.Open(infile, 1) # 0 means read-only. 1 means writeable.\n",
  105. " layer = dataSource.GetLayer()\n",
  106. " sourceprj = layer.GetSpatialRef()\n",
  107. " transform = osr.CoordinateTransformation(sourceprj, targetprj)\n",
  108. " \n",
  109. " # Create the output shapefile\n",
  110. " outDriver = ogr.GetDriverByName(\"Esri Shapefile\")\n",
  111. " outDataSource = outDriver.CreateDataSource(outfile)\n",
  112. " outlayer = outDataSource.CreateLayer('', targetprj, ogr.wkbPolygon)\n",
  113. " outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))\n",
  114. " \n",
  115. " #Iterate over Features\n",
  116. " i = 0\n",
  117. " for feature in layer:\n",
  118. " transformed = feature.GetGeometryRef()\n",
  119. " transformed.Transform(transform) #reproject geometry\n",
  120. "\n",
  121. " geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb()) # create geometry from wkb (write geometry of reprojected geometry)\n",
  122. " defn = outlayer.GetLayerDefn() #layer definition\n",
  123. " feat = ogr.Feature(defn) #create new feature\n",
  124. " feat.SetField('id', i) #set id\n",
  125. " feat.SetGeometry(geom) #set geometry\n",
  126. " outlayer.CreateFeature(feat) \n",
  127. " i += 1\n",
  128. " feat = None\n",
  129. " \n",
  130. "def array2raster(array, geoTransform, projection, filename):\n",
  131. " \"\"\" \n",
  132. " This function tarnsforms a numpy array to a geotiff projected raster\n",
  133. " input: array array (n x m) input array\n",
  134. " geoTransform tuple affine transformation coefficients\n",
  135. " projection string projection\n",
  136. " filename string output filename\n",
  137. " output: dataset gdal raster dataset\n",
  138. " dataset.GetRasterBand(1) band object of dataset\n",
  139. " \n",
  140. " \"\"\"\n",
  141. " pixels_x = array.shape[1]\n",
  142. " pixels_y = array.shape[0]\n",
  143. " \n",
  144. " driver = gdal.GetDriverByName('GTiff')\n",
  145. " dataset = driver.Create(\n",
  146. " filename,\n",
  147. " pixels_x,\n",
  148. " pixels_y,\n",
  149. " 1,\n",
  150. " gdal.GDT_Float64, )\n",
  151. " dataset.SetGeoTransform(geoTransform)\n",
  152. " dataset.SetProjection(projection)\n",
  153. " dataset.GetRasterBand(1).WriteArray(array)\n",
  154. " dataset.FlushCache() # Write to disk.\n",
  155. " return dataset, dataset.GetRasterBand(1) #If you need to return, remenber to return also the dataset because the band don`t live without dataset.\n",
  156. " \n",
  157. "def clip_raster(filename, shp):\n",
  158. " \"\"\"\n",
  159. " This function clips a raster based on a shapefile\n",
  160. " input: filename string input raster filename\n",
  161. " shp dataframe input shapefile open with geopandas\n",
  162. " output: clipped array (1 x n x m) clipped array \n",
  163. " clipped_meta dict metadata\n",
  164. " cr_ext tuple extent of clipped data\n",
  165. " gt tuple affine transformation coefficients\n",
  166. " \"\"\"\n",
  167. " inraster = rasterio.open(filename)\n",
  168. " \n",
  169. " extent_geojson = mapping(shp['geometry'][0])\n",
  170. " clipped, crop_affine = mask(inraster, \n",
  171. " shapes=[extent_geojson], \n",
  172. " nodata = np.nan,\n",
  173. " crop=True)\n",
  174. " clipped_meta = inraster.meta.copy()\n",
  175. " # Update the metadata to have the new shape (x and y and affine information)\n",
  176. " clipped_meta.update({\"driver\": \"GTiff\",\n",
  177. " \"height\": clipped.shape[0],\n",
  178. " \"width\": clipped.shape[1],\n",
  179. " \"transform\": crop_affine})\n",
  180. " cr_ext = rasterio.transform.array_bounds(clipped_meta['height'], \n",
  181. " clipped_meta['width'], \n",
  182. " clipped_meta['transform'])\n",
  183. " \n",
  184. " # transform to gdal\n",
  185. " gt = crop_affine.to_gdal()\n",
  186. " \n",
  187. " return clipped, clipped_meta, cr_ext, gt\n",
  188. " \n",
  189. "def reclassify(array):\n",
  190. " \"\"\"\n",
  191. " This function reclassifies an array\n",
  192. " input: array array (n x m) input array\n",
  193. " output: reclass array (n x m) reclassified array\n",
  194. " \"\"\"\n",
  195. " reclass = np.zeros((array.shape[0],array.shape[1]))\n",
  196. " for i in range(0,array.shape[0]):\n",
  197. " for j in range(0,array.shape[1]):\n",
  198. " if math.isnan(array[i,j]):\n",
  199. " reclass[i,j] = np.nan\n",
  200. " elif array[i,j] < 0.1:\n",
  201. " reclass[i,j] = 1\n",
  202. " elif array[i,j] < 0.27:\n",
  203. " reclass[i,j] = 2\n",
  204. " elif array[i,j] < 0.44:\n",
  205. " reclass[i,j] = 3\n",
  206. " elif array[i,j] < 0.66:\n",
  207. " reclass[i,j] = 4\n",
  208. " else:\n",
  209. " reclass[i,j] = 5\n",
  210. " \n",
  211. " return reclass\n",
  212. " "
  213. ]
  214. },
  215. {
  216. "cell_type": "markdown",
  217. "metadata": {},
  218. "source": [
  219. "In the following cell you should have to modify the paths for the images and the shapefile"
  220. ]
  221. },
  222. {
  223. "cell_type": "code",
  224. "execution_count": null,
  225. "metadata": {},
  226. "outputs": [],
  227. "source": [
  228. "# Paths\n",
  229. "path_prefire = \"F:/Burn_Severity/S2A_MSIL2A_20161220T143742_N0204_R096_T18HYF_20161220T145131.SAFE/GRANULE/L2A_T18HYF_A007815_20161220T145131/IMG_DATA/R20m/\"\n",
  230. "path_postfire = \"F:/Burn_Severity/S2A_MSIL2A_20170218T143751_N0204_R096_T18HYF_20170218T145150.SAFE/GRANULE/L2A_T18HYF_A008673_20170218T145150/IMG_DATA/R20m/\"\n",
  231. "# Define shapefile\n",
  232. "infile_shp = \"F:/Burn_Severity/Empedrado_adm_boundary/Empedrado.shp\"\n",
  233. "# Define reprojected shapefile\n",
  234. "outfile_shp = \"F:/Burn_Severity/Empedrado_adm_boundary/projected.shp\"\n",
  235. "# name of the output dNBR raster file\n",
  236. "filename = \"F:/Burn_Severity/dNBR.tiff\"\n",
  237. "# name of clipped dNBR raster\n",
  238. "filename2 = \"F:/Burn_Severity/dNBR_clipped.tiff\"\n",
  239. "# path to save figure\n",
  240. "fname = \"F:/Burn_Severity/map.png\"\n",
  241. "\n",
  242. "# Sentinel-2 Bands used for NBR calculation \n",
  243. "band1 = '8A'\n",
  244. "band2 = '12'\n",
  245. " \n",
  246. "# Read the pre-fire band images \n",
  247. "(pre_fire_b8a, crs, geoTransform, targetprj) = read_band_image(band1, path_prefire)\n",
  248. "(pre_fire_b12, crs, geoTransform, targetprj) = read_band_image(band2, path_prefire)\n",
  249. " \n",
  250. "# Calculation of pre-fire NBR\n",
  251. "pre_fire_nbr = nbr(pre_fire_b8a.astype(int),pre_fire_b12.astype(int))\n",
  252. "\n",
  253. "# Read the post-fire band images\n",
  254. "(post_fire_b8a, crs, geoTransform, targetprj) = read_band_image(band1, path_postfire)\n",
  255. "(post_fire_b12, crs, geoTransform, targetprj) = read_band_image(band2, path_postfire)\n",
  256. " \n",
  257. "# Calculation of post-fire NBR\n",
  258. "post_fire_nbr = nbr(post_fire_b8a.astype(int),post_fire_b12.astype(int))\n",
  259. " \n",
  260. "# Calculation of dNBR\n",
  261. "DNBR = dnbr(pre_fire_nbr,post_fire_nbr)\n",
  262. " \n",
  263. "# Reprojection of shapefile with gdal to match projection of Sentinel-2 images\n",
  264. "reproject_shp_gdal(infile_shp, outfile_shp, targetprj)\n",
  265. " \n",
  266. "# Read the reprojected shapefile\n",
  267. "fire_boundary = gpd.read_file(outfile_shp)\n",
  268. " \n",
  269. "# project dNBR to images projection\n",
  270. "dnbr_tif, dnbr_tifBand = array2raster(DNBR, geoTransform, crs, filename)\n",
  271. " \n",
  272. "# clip raster dNBR file to Empedrado shapefile\n",
  273. "(clipped_dnbr, clipped_dnbr_meta, cr_extent, gt) = clip_raster(filename, fire_boundary)\n",
  274. "clipped_ds , clipped_ds_rasterband = array2raster(clipped_dnbr[0], gt, crs, filename2)\n"
  275. ]
  276. },
  277. {
  278. "cell_type": "markdown",
  279. "metadata": {},
  280. "source": [
  281. "Plot the results"
  282. ]
  283. },
  284. {
  285. "cell_type": "code",
  286. "execution_count": null,
  287. "metadata": {},
  288. "outputs": [],
  289. "source": [
  290. "# plot \n",
  291. "\n",
  292. "#set colors for plotting and classes\n",
  293. "cmap = matplotlib.colors.ListedColormap(['green','yellow','orange','red','purple'])\n",
  294. "cmap.set_over('purple')\n",
  295. "cmap.set_under('white')\n",
  296. "bounds = [-0.5, 0.1, 0.27, 0.440, 0.660, 1.3] \n",
  297. "norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N) \n",
  298. "\n",
  299. "fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'xticks': [], 'yticks': []})\n",
  300. "cax = ax.imshow(clipped_ds_rasterband.ReadAsArray(), cmap=cmap, norm = norm)\n",
  301. "plt.title('Burn Severity Map')\n",
  302. "cbar = fig.colorbar(cax, ax=ax, fraction=0.035, pad=0.04, ticks=[-0.2, 0.18, 0.35, 0.53, 1])\n",
  303. "cbar.ax.set_yticklabels(['Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity'])\n",
  304. "plt.show()\n",
  305. "plt.savefig(fname, bbox_inches=\"tight\") \n",
  306. " "
  307. ]
  308. },
  309. {
  310. "cell_type": "markdown",
  311. "metadata": {},
  312. "source": [
  313. "Calculate burnt area"
  314. ]
  315. },
  316. {
  317. "cell_type": "code",
  318. "execution_count": null,
  319. "metadata": {},
  320. "outputs": [],
  321. "source": [
  322. "# calculate burnt area (pixel size 20m*20m)\n",
  323. "reclass = reclassify(clipped_ds_rasterband.ReadAsArray())\n",
  324. "k = ['Unburned hectares', 'Low severity hectares', 'Moderate-low severity hectares', 'Moderate-high severity hectares', 'High severity']\n",
  325. "for i in range(1,6):\n",
  326. " x = reclass[reclass == i]\n",
  327. " l= x.size*0.04\n",
  328. " print(\"%s: %.2f\" % (k[i-1], l))"
  329. ]
  330. },
  331. {
  332. "cell_type": "code",
  333. "execution_count": null,
  334. "metadata": {},
  335. "outputs": [],
  336. "source": []
  337. }
  338. ],
  339. "metadata": {
  340. "kernelspec": {
  341. "display_name": "Python 3",
  342. "language": "python",
  343. "name": "python3"
  344. },
  345. "language_info": {
  346. "codemirror_mode": {
  347. "name": "ipython",
  348. "version": 3
  349. },
  350. "file_extension": ".py",
  351. "mimetype": "text/x-python",
  352. "name": "python",
  353. "nbconvert_exporter": "python",
  354. "pygments_lexer": "ipython3",
  355. "version": "3.6.5"
  356. }
  357. },
  358. "nbformat": 4,
  359. "nbformat_minor": 2
  360. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!