|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +""" firemap |
| 3 | +===== |
| 4 | +This package provides functions to download annual fire frequencies from |
| 5 | +Google Earth Engine (i.e. number of months with active fires), based on |
| 6 | +monthly, global burned area maps derived with MODIS at a 500 m resolution. |
| 7 | +For details on the dataset, refer to the GEE's data catalog |
| 8 | +(https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD64A1). |
| 9 | +The download output is then used to derives spatial-temporal metrics that |
| 10 | +characterize pixel-level fire regimes: |
| 11 | +
|
| 12 | +The algorithm evaluates the fire regime of an area given GeoTiff images on |
| 13 | +annual fire occurrences. For every pixel, the algorithm will estimate the |
| 14 | +following metrics: |
| 15 | + |
| 16 | + * Fire Return Interval (FRI) - Mean number of years between firest |
| 17 | + * Min. return time (MRI1) - Minimum between-burn interval |
| 18 | + * Mean return time (MRI2) - Mean between-burn interval |
| 19 | + * Maximum return time (MRI3) - Maximum between-burn interval |
| 20 | +
|
| 21 | +The FRI is expressed as the quotient between the number of years with fires |
| 22 | +and the number of years in the time-series. Then, for each pixel, fire map |
| 23 | +uses Running Length Encoding (RLE) to break a time-series of fire occurrences |
| 24 | +(with 1 for "fire" and 0 for "unburnt") into segments of equal value. The |
| 25 | +length of segments corresponding to "unburnt" periods are used to calculate |
| 26 | +MRI1, MRI2, and MRI3. |
| 27 | +Note: the algorithm was designed to download data with a maximum resolution |
| 28 | +of 1-km, which is limited by the quota imposed by GEE on the size of files |
| 29 | +that can be directly downloaded after processing. Applying this algorithm to |
| 30 | +larger files will require some modifications to the gee_download() function. |
| 31 | +""" |
| 32 | + |
| 33 | +from progress.bar import Bar |
| 34 | +import urllib.request as ur |
| 35 | +from os.path import join |
| 36 | +from osgeo import gdal |
| 37 | +import rasterio as rt |
| 38 | +import numpy as np |
| 39 | +import glob2 as g |
| 40 | +import shutil |
| 41 | +import rle |
| 42 | +import os |
| 43 | +import ee |
| 44 | + |
| 45 | +def gee_download(data_path, year): |
| 46 | + |
| 47 | + """ Returns annaul |
| 48 | + |
| 49 | + Sum and return the ages of your son and daughter |
| 50 | + |
| 51 | + Parameters |
| 52 | + ------------ |
| 53 | + path: str |
| 54 | + Path to directory where to download the data |
| 55 | + path: str |
| 56 | + Year for which to download data |
| 57 | + Return |
| 58 | + ----------- |
| 59 | + age : int |
| 60 | + The sum of your son and daughter ages. |
| 61 | + """ |
| 62 | + |
| 63 | + # connect to GEE |
| 64 | + ee.Initialize() |
| 65 | + |
| 66 | + #------------------------------------------------------------------------# |
| 67 | + bar = Bar('downloading data from GEE', max=162) |
| 68 | + #------------------------------------------------------------------------# |
| 69 | + |
| 70 | + ydir = os.path.join(data_path, year) |
| 71 | + if os.path.isdir(ydir) == False: |
| 72 | + os.mkdir(ydir) |
| 73 | + |
| 74 | + for x in range(0,18): |
| 75 | + |
| 76 | + for y in range(0,9): |
| 77 | + |
| 78 | + #================================================================# |
| 79 | + # |
| 80 | + #================================================================# |
| 81 | + |
| 82 | + # build name of zip file to download |
| 83 | + zfile = ydir + "{:02d}".format(x) + '-' + \ |
| 84 | + "{:02d}".format(y) + '_' + year + '.zip' |
| 85 | + |
| 86 | + # build name of zip file to download |
| 87 | + zdir = ydir + "{:02d}".format(x) + '-' + \ |
| 88 | + "{:02d}".format(y) + '_' + year |
| 89 | + |
| 90 | + #================================================================# |
| 91 | + # |
| 92 | + #================================================================# |
| 93 | + |
| 94 | + # output extent |
| 95 | + xmin = -180+(20*x) |
| 96 | + ymin = 90-(20*(y+1)) |
| 97 | + xmax= -180+(20*(x+1)) |
| 98 | + ymax = 90-(20*y) |
| 99 | + |
| 100 | + # boulding box to write |
| 101 | + geometry = ('[[' + str(xmin) + ' ,' + str(ymax) + |
| 102 | + '], [' + str(xmax) + ', ' + str(ymax) + |
| 103 | + '], [' + str(xmax) + ', ' + str(ymin) + |
| 104 | + '], [' + str(xmin) + ', ' + str(ymin) + ']]') |
| 105 | + |
| 106 | + # subset geometry |
| 107 | + geometry = ee.Geometry.Rectangle(xmin,ymin,xmax,ymax) |
| 108 | + |
| 109 | + #================================================================# |
| 110 | + # |
| 111 | + #================================================================# |
| 112 | + |
| 113 | + # read collection |
| 114 | + col = ee.ImageCollection('MODIS/006/MCD64A1').filterDate( |
| 115 | + year + '-01-01', year + '-12-31' |
| 116 | + ).filterBounds(geometry).select('BurnDate') |
| 117 | + |
| 118 | + #================================================================# |
| 119 | + # |
| 120 | + #================================================================# |
| 121 | + |
| 122 | + # boulding box to write |
| 123 | + geometry = ('[[' + str(xmin) + ' ,' + str(ymax) + |
| 124 | + '], [' + str(xmax) + ', ' + str(ymax) + |
| 125 | + '], [' + str(xmax) + ', ' + str(ymin) + |
| 126 | + '], [' + str(xmin) + ', ' + str(ymin) + ']]') |
| 127 | + |
| 128 | + def main(image): |
| 129 | + return image.expression('b("BurnDate") > 0').rename('burn') |
| 130 | + |
| 131 | + image = col.map(main).sum() |
| 132 | + |
| 133 | + # extract download link |
| 134 | + path = image.getDownloadUrl({'description':year, 'name':year, |
| 135 | + 'reducer':'mean', 'scale':1000, |
| 136 | + 'crs':'EPSG:4326', 'region':geometry}) |
| 137 | + |
| 138 | + #================================================================# |
| 139 | + # |
| 140 | + #================================================================# |
| 141 | + |
| 142 | + # download file |
| 143 | + while os.path.isfile(zfile) == False: |
| 144 | + try: |
| 145 | + ur.urlretrieve(path, zfile) |
| 146 | + except: |
| 147 | + continue |
| 148 | + |
| 149 | + # create unziped directory |
| 150 | + if os.path.isfile(zdir) == False: |
| 151 | + os.mkdir(zdir) |
| 152 | + |
| 153 | + # unpack file into new directory |
| 154 | + shutil.unpack_archive(zfile, zdir) |
| 155 | + |
| 156 | + bar.next() |
| 157 | + print(str(x) + '-' + str(y)) |
| 158 | + |
| 159 | + bar.next() |
| 160 | + |
| 161 | + bar.finish() |
| 162 | + |
| 163 | + #------------------------------------------------------------------------# |
| 164 | + bar = Bar('mosaic downloaded files', max=162) |
| 165 | + #------------------------------------------------------------------------# |
| 166 | + |
| 167 | + # list files |
| 168 | + files = g.glob(join(data_path, '**', '*.tif')) |
| 169 | + bar.next() |
| 170 | + |
| 171 | + iname = join(data_path, year, '.vrt') |
| 172 | + gdal.BuildVRT(iname, files) |
| 173 | + bar.next() |
| 174 | + |
| 175 | + oname = join(data_path, year, '.tif') |
| 176 | + gdal.Translate(iname, oname) |
| 177 | + bar.next() |
| 178 | + bar.finish() |
| 179 | + |
| 180 | + |
| 181 | + |
| 182 | +def fire_regime(input_path, output_path): |
| 183 | + |
| 184 | + #------------------------------------------------------------------------# |
| 185 | + # list input images |
| 186 | + #------------------------------------------------------------------------# |
| 187 | + |
| 188 | + files = sorted(g.glob(join(input_path, '*.tif'))) |
| 189 | + |
| 190 | + #------------------------------------------------------------------------# |
| 191 | + # use first image as a reference and extract metadata |
| 192 | + #------------------------------------------------------------------------# |
| 193 | + |
| 194 | + p = rt.open(files[0]).meta.copy() |
| 195 | + p.update(dtype='float32', compress='deflate', predict=2, zlevel=9) |
| 196 | + |
| 197 | + #------------------------------------------------------------------------# |
| 198 | + bar = Bar('estimate Fire Recurrence Interval (FRI)', max=len(files)+2) |
| 199 | + #------------------------------------------------------------------------# |
| 200 | + |
| 201 | + # estimate sum of years with fires |
| 202 | + ia = np.zeros((p['height'],p['width']), dtype='uint8') |
| 203 | + for f in files: |
| 204 | + ia = ia + (rt.open(f).read(1) > 0).astype('uint8') |
| 205 | + bar.next() |
| 206 | + |
| 207 | + # find non-zero pixels |
| 208 | + px = np.where(ia > 0) |
| 209 | + |
| 210 | + # estimate average interval |
| 211 | + oa = np.zeros(ia.shape, dtype='float32') |
| 212 | + oa[px] = ia[px] / len(files) |
| 213 | + bar.next() |
| 214 | + |
| 215 | + #------------------------------------------------------------------------# |
| 216 | + # write fire regime image |
| 217 | + #------------------------------------------------------------------------# |
| 218 | + |
| 219 | + ods = rt.open(join(output_path, 'fire_recurrence_interval.tif'), 'w', **p) |
| 220 | + ods.write(oa, indexes=1) |
| 221 | + ods.close() |
| 222 | + bar.next() |
| 223 | + bar.finish() |
| 224 | + |
| 225 | + #------------------------------------------------------------------------# |
| 226 | + bar = Bar('evaluate intervals between fires', max=len(px[0])+len(files)+3) |
| 227 | + #------------------------------------------------------------------------# |
| 228 | + |
| 229 | + # use non-zero pixel indices to derive a new matrix with all years |
| 230 | + ia = np.zeros((len(px[0]),len(files))) |
| 231 | + for i in range(0,len(files)): |
| 232 | + ia[:,i] = rt.open(files[i]).read(1)[px] |
| 233 | + bar.next() |
| 234 | + |
| 235 | + # create empty ouput array |
| 236 | + oa = np.zeros((len(px[0]),3)) |
| 237 | + |
| 238 | + # derive metrics on the longest time without fires |
| 239 | + for i in range(0,len(px[0])): |
| 240 | + unique_id, id_length = rle.encode(ia[i,:]) # find sequences |
| 241 | + ind = np.where(np.array(unique_id) == 0) |
| 242 | + if len(ind[0]) > 0: |
| 243 | + oa[i,0] = np.min(np.array(id_length)[ind]) # min. interval |
| 244 | + oa[i,1] = np.max(np.array(id_length)[ind]) # mean interval |
| 245 | + oa[i,2] = np.mean(np.array(id_length)[ind]) # max. interval |
| 246 | + |
| 247 | + bar.next() |
| 248 | + |
| 249 | + #------------------------------------------------------------------------# |
| 250 | + # write metrics |
| 251 | + #------------------------------------------------------------------------# |
| 252 | + |
| 253 | + # update metadata to reduce file size |
| 254 | + p.update(dtype='uint8') |
| 255 | + |
| 256 | + # assign estimate metric to output array |
| 257 | + fr = np.zeros((p['height'],p['width']), dtype='uint8') |
| 258 | + fr[px] = oa[:,0] |
| 259 | + |
| 260 | + # write output |
| 261 | + ods = rt.open(join(input_path, 'mininum_return.tif'), 'w', **p) |
| 262 | + ods.write(fr, indexes=1) |
| 263 | + ods.close() |
| 264 | + bar.next() |
| 265 | + |
| 266 | + # assign estimate metric to output array |
| 267 | + fr = np.zeros((p['height'],p['width']), dtype='uint8') |
| 268 | + fr[px] = oa[:,1] |
| 269 | + |
| 270 | + # write output |
| 271 | + ods = rt.open(join(input_path, 'maximum_return.tif'), 'w', **p) |
| 272 | + ods.write(fr, indexes=1) |
| 273 | + ods.close() |
| 274 | + bar.next() |
| 275 | + |
| 276 | + # assign estimate metric to output array |
| 277 | + fr = np.zeros((p['height'],p['width']), dtype='uint8') |
| 278 | + fr[px] = oa[:,2] |
| 279 | + |
| 280 | + # write output |
| 281 | + ods = rt.open(join(input_path, 'mean_return.tif'), 'w', **p) |
| 282 | + ods.write(fr, indexes=1) |
| 283 | + ods.close() |
| 284 | + bar.next() |
| 285 | + bar.finish() |
0 commit comments