Fire Map Stacking

Fire Map Stacking.txt

#Fire Map Stacking
#Version 1.0
#Xinyuan Wei 12/20/2015
#Fire Cycle Study
'''
This program can stack simulated fire maps and create a TSLF map.
'''
import os.path
import glob
import gdal
import osr
import ogr
import numpy as np
from gdalconst import *
from PIL import Image
gdal.AllRegister()

#create the base map (array)
dataname='1.img'
rasterdata=gdal.Open(dataname,GA_ReadOnly)

cols=rasterdata.RasterXSize
rows=rasterdata.RasterYSize
bands=rasterdata.RasterCount
driver=rasterdata.GetDriver().LongName
    
geotransform=rasterdata.GetGeoTransform()
originX=geotransform[0]
originY=geotransform[3]
north=geotransform[4]
pixelWidth=geotransform[1]
pixelHeight=geotransform[5]

#the function converts an array to a raster map
def arraytoraster(newRaster,rasterOrigin,pixelWidth,pixelHeight,array):
    cols=array.shape[1]
    rows=array.shape[0]
    originX=rasterOrigin[0]
    originY=rasterOrigin[1]

    driver=gdal.GetDriverByName('GTiff')
    outRaster=driver.Create(newRaster,cols,rows,1,gdal.GDT_UInt16)
    outRaster.SetGeoTransform((originX,pixelWidth,0,originY,0,pixelHeight))
    outband=outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS=osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

#set the initial TSLF in each grid as 1200
band=rasterdata.GetRasterBand(1)  
base_map=band.ReadAsArray(0,0,cols,rows)
for i in range (cols):
    for j in range(rows):
        base_map[i,j]=1200

n=0 #create test maps
#stack fire maps
for image in glob.glob('*.img'):
    print(image)
    fire_map=gdal.Open(image,GA_ReadOnly)
    band_map=fire_map.GetRasterBand(1)    
    fire_year=band_map.ReadAsArray(0,0,cols,rows)

    for i in range (cols):
        for j in range(rows):
            if fire_year[i,j]!=0:
                if base_map[i,j]>fire_year[i,j]:
                    base_map[i,j]=fire_year[i,j]

    n=n+1
    print(n)
    if n%20==0:
        rasterOrigin=(originX,originY)
        newRaster='TSLF_%d.img'%n
        array=np.array(base_map)
        
        arraytoraster(newRaster,rasterOrigin,pixelWidth,pixelHeight,array)