#!/usr/bin/env python import os import sys import h5py import argparse import numpy as np import datetime from PIL import Image from pyunpack import Archive from rich.progress import track from rich import print """geraFocos.py: Como usar: faça ./geraFocos.py ex: ./geraFocos.py """ __author__ = "Luiz Flávio Rodrigues" __copyright__ = "Copyright 2020, INPE" __credits__ = ["Luiz Flavio Rodrigues"] __license__ = "CC-GPL" __version__ = "1.0" __maintainer__ = "Luiz Flávio Rodrigues" __email__ = "luiz.rodrigues@inpe.br" __status__ = "Production" def main(args): ano=int(args[1]) mes=int(args[2]) dia=int(args[3]) viirs = args[4] ecmwf = args[5] inpe = args[6] meteosat = args[7] print("") print("[bold green] ____ ____ _ __ __ ____ [/bold green]") print("[bold green]| __ )| _ \ / \ | \/ / ___| [/bold green]") print("[bold green]| _ \| |_) | / _ \ | |\/| \___ \ [/bold green]") print("[bold green]| |_) | _ < / ___ \| | | |___) |[/bold green]") print("[bold green]|____/|_| \_\/_/ \_\_| |_|____/ [/bold green]") print("[bold green] PRE - PREP - Gera de Focos de Fogo[/bold green]") print("Rev 1.0 - e-mail: luiz.rodrigues@inpe.br") print("Author: Luiz Flávio Rodrigues INPE:copyright:") """ ##################################################### ____ ____ _ | _ \/ ___| / \ | | | \___ \ / _ \ | |_| |___) / ___ \ |____/|____/_/ \_\ ###################################################### """ if inpe=="s" or inpe=="S": print("") print("Processando dados DSA:") dsa_input_dir = "/share/bramsrd/dist/BRAMS/data/EMISSIONS/DSA/" dsa_output_dir = "/share/bramsrd/dist/BRAMS/data/FOCOS/{0:04d}{1:02d}{2:02d}/DSA/".format(ano,mes,dia) dsa_file_name = dsa_input_dir+"Focos{0:04d}{1:02d}{2:02d}.txt".format(ano,mes,dia) dsa_focos_file_name = dsa_output_dir+"Focos.{0:04d}{1:02d}{2:02d}.proc.txt".format(ano,mes,dia) dsa_file = open(dsa_file_name,"r") dsa_focos_file = open(dsa_focos_file_name,"w") valid_sensor = ["NOAA-18D","NOAA-19D","NOAA-18","NOAA-19","NOAA-15D","NOAA-15","NPP_375","METOP-B","TERRA_M-M","TERRA_M-T","AQUA_M-M","AQUA_M-T"] lines = dsa_file.readlines() header = True count = 0 for line in track(lines): if header: header=False continue fields = line.split(",") lat = float(fields[2]) lon = float(fields[3]) sensor = fields[6] municipio = fields[7] estado = fields[8] pais = fields[9] if sensor in valid_sensor: if lon>-100.0 and lon<-30.: if lat>-60.0 and lat<15: count = count+1 dsa_focos_file.write("{0:12.8f} {1:12.8f} ".format(lat,lon)+municipio+"-"+estado+"-"+pais+"\n") dsa_file.close() dsa_focos_file.close() print("Processados {0:6d} pontos para dados da DSA!".format(count)) print("Escritos no arquivo "+dsa_focos_file_name) """ ###################################### __ _____ ___ ____ ____ \ \ / /_ _|_ _| _ \/ ___| \ \ / / | | | || |_) \___ \ \ V / | | | || _ < ___) | \_/ |___|___|_| \_\____/ ###################################### """ if viirs=="s" or viirs=="S": print("") print("Processando dados VIIRS:") viirs_input_dir = "/share/bramsrd/dist/BRAMS/data/EMISSIONS/VIIRS/{0:04d}{1:02d}{2:02d}/".format(ano,mes,dia) viirs_output_dir = "/share/bramsrd/dist/BRAMS/data/FOCOS/{0:04d}{1:02d}{2:02d}/MODIS/".format(ano,mes,dia) viirs_file_name = os.listdir(viirs_input_dir) viirs_focos_file_name = viirs_output_dir+"Fires.{0:04d}{1:02d}{2:02d}.txt".format(ano,mes,dia) viirs_focos_file = open(viirs_focos_file_name,"w") count = 0 for file in track(viirs_file_name): viirs_file = open(viirs_input_dir+file,"r") lines = viirs_file.readlines() header = True for line in lines: if header: header = False continue fields = line.split(",") lat = float(fields[5]) lon = float(fields[6]) frp = float(fields[10]) if frp>0.0: if lon>-100.0 and lon<-30.: if lat>-60.0 and lat<15: count = count+1 viirs_focos_file.write("{0:12.8f} {1:12.8f} \n ".format(lon,lat)) viirs_file.close() viirs_focos_file.close() print("Processados {0:6d} pontos para dados viirs!".format(count)) print("Escritos no arquivo "+viirs_focos_file_name) """ __ __ _____ _____ _____ ___ ____ _ _____ | \/ | ____|_ _| ____/ _ \/ ___| / \|_ _| | |\/| | _| | | | _|| | | \___ \ / _ \ | | | | | | |___ | | | |__| |_| |___) / ___ \| | |_| |_|_____| |_| |_____\___/|____/_/ \_\_| """ if meteosat == "s" or meteosat == "S": print("") print("Processando dados Meteosat:") meteosat_input_dir = "/share/bramsrd/dist/BRAMS/data/EMISSIONS/METEOSAT/{0:04d}/{1:02d}/".format(ano,mes) meteosat_output_dir = "/share/bramsrd/dist/BRAMS/data/FOCOS/{0:04d}{1:02d}{2:02d}/EXTRA/".format(ano,mes,dia) meteosat_focos_file_name = meteosat_output_dir+"Fires.extra.{0:04d}{1:02d}{2:02d}.txt".format(ano,mes,dia) meteosat_focos_file = open(meteosat_focos_file_name,"w") count=0 #meteosat_file_name = os.listdir(meteosat_input_dir+"HDF5_LSASAF_MSG_FRP-PIXEL-ListProduct_MSG-Disk_{0:04d}{1:02d}{2:02d}*".format(ano,mes,dia)) file_list=[] for h in range(24): for m in range(0,60,15): file_list.append(meteosat_input_dir+"HDF5_LSASAF_MSG_FRP-PIXEL-ListProduct_MSG-Disk_{0:04d}{1:02d}{2:02d}{3:02d}{4:02d}".format(ano,mes,dia,h,m)) for meteosat_file_name in track(file_list): present = os.path.exists(meteosat_file_name) if not present: if not os.path.exists(meteosat_file_name+".bz2"): print("Arquivo "+meteosat_file_name+" não existe nos formatos hdf5 ou bzip2!") continue else: Archive(meteosat_file_name+".bz2").extractall(meteosat_input_dir) try: f = h5py.File(meteosat_file_name,'r') except: print("problema: %s não parece ser um arquivo hdf5 válido" % meteosat_file_name) n1 = f.get('FRP') n1 = np.array(n1) frp = n1.tolist() n2 = f.get('LATITUDE') print(n2) n2 = np.array(n2)/100. print(n2) lat = n2.tolist() print(lat) n3 = f.get('LONGITUDE') n3 = np.array(n3)/100. lon = n3.tolist() n4 = f.get('PIXEL_SIZE') n4 = np.array(n4)/100. psiz = n4.tolist() n5 = f.get('FIRE_CONFIDENCE') n5 = np.array(n5)/100. conf = n5.tolist() for pos in range(len(frp)): if conf[pos]>0.85: if frp[pos]>1.0: if lon[pos]>-90 and lon[pos]<-30.0: if lat[pos]>-60.0 and lat[pos]<15: meteosat_focos_file.write("{0:12.8f} {1:12.8f} {2:12.8f} \n ".format(lon[pos],lat[pos],psiz[pos])) count=count+1 f.close() meteosat_focos_file.close() print("Processados {0:6d} pontos para dados meteosat!".format(count)) print("Escritos no arquivo "+meteosat_focos_file_name) """ ____ _ __ __ ____ ____ ___ _____ ____ / ___| / \ | \/ / ___| / ___|/ _ \| ____/ ___| | | / _ \ | |\/| \___ \ _____| | _| | | | _| \___ \ | |___ / ___ \ | | | |___) |_____| |_| | |_| | |___ ___) | \____/_/ \_\_| |_|____/ \____|\___/|_____|____/ """ if ecmwf=="s" or ecmwf == "S": print("") #Le o arquivo de vegetação do USGS print("Pegando tiff de vegetação USGS...") Image.MAX_IMAGE_PIXELS = None im = Image.open('/share/bramsrd/dist/BRAMS/data/GEOFIRE/gbogegeo20.tif') imarray = np.array(im) nLats=imarray.shape[0] nLons=imarray.shape[1] print("Pegando tamanhos do TIFF ...") fSizes=open("/share/bramsrd/dist/BRAMS/data/GEOFIRE/gbogegeo20.tfw","r") info=fSizes.readlines() #Obtém o deltaLat e deltaLon do arquivo da USGS delta=float(info[0]) #Cria uma lista de fire Área com 100 elementos (1 para cada bioma) #E preenche com valor default 0.050 km2 fireArea=[0.050]*200 #Le a climatologia de tamanho de fogo por bioma print("Pegando dados de climatologia por ECO ...") fEco=open("/share/bramsrd/dist/BRAMS/data/GEOFIRE/fireClimatology.txt","r") ecoInfo=fEco.readlines() #Obtem o valor da fire área por bioma da climatologia for line in ecoInfo: fireArea[int(line.split(";")[0])]=float(line.split(";")[3])/1.0e+6 #Determina as latitudes e longitudes de cada ponto do arquivo USGS lat = [] lon = [] for i in range(nLats): lat.append(90.0 - i * delta) for j in range(nLons): lon.append(-180 + j * delta) print("") print("Processando dados GOES-CAMS:") horas=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] minutos=[0,10,20,30,40,50] fileErr="/share/bramsrd/dist/BRAMS/data/FOCOS/{0:04d}{1:02d}{2:02d}/GOES-CAMS/FileErr.txt".format(ano,mes,dia) ferr=open(fileErr,"w") file_list=[] for h in horas: for m in minutos: file_list.append("/share/bramsrd/dist/BRAMS/data/EMISSIONS/GOES-ECMWF/{0:04d}{1:02d}{2:02d}/CAMS__GOES16_FRP-PIXEL-ListProduct_GOESE-FD_{0:04d}{1:02d}{2:02d}{3:02d}{4:02d}.h5".format(ano,mes,dia,h,m)) npoints=0 for fileIn in track(file_list): x=fileIn.split('_') pos=len(x)-1 #print("Pos {0:2d}".format(pos)) #print("X[pos]= %s" % x[pos]) x=x[pos].split('.') data=x[0] #print(data[8:10]) #print(data[10:12]) hora=int(data[8:10])#h min=int(data[10:12])#m #print("Hora: {0:2d}, MIn: {1:2d}".format(hora,min)) a=datetime.datetime.strptime(data,'%Y%m%d%H%M') djul=int(a.strftime('%j')) #check if the file passed as agument exists if (not os.path.isfile(fileIn)): print("problem: the file %s doesnt exist." % fileIn) ferr.write("problem: the file %s doesnt exist." % fileIn) continue try: #open the file in reading mode f = h5py.File(fileIn,'r') except: print("problem: %s is not a valid hdf5 file" % fileIn) ferr.write("problem: %s is not a valid hdf5 file" % fileIn) continue n1 = f.get('FRP') n1=np.array(n1) if n1.any()==None: ferr.write(fileIn) continue n1=n1/10.0 n1=n1.tolist() #print(n1) n2 = f.get('LATITUDE') escala=n2.attrs["SCALING_FACTOR"] n2=np.array(n2)/escala n2=n2.tolist() #print(n2) n3 = f.get('LONGITUDE') escala=n3.attrs["SCALING_FACTOR"] n3=np.array(n3)/escala n3=n3.tolist() # n4 = f.get('PIXEL_VZA') n4=np.array(n4)/100. n4=n4.tolist() #print(n3) n5 = f.get('PIXEL_SIZE') n5=np.array(n5)/100. n5=n5*1000 #Para metros n5=n5.tolist() #print(n3) n6 = f.get('BT_MIR') n6=np.array(n6)/10. n6=n6.tolist() #print(n3) n7 = f.get('BT_TIR') n7=np.array(n7)/10. n7=n7.tolist() #print(n3) n8 = f.get('FIRE_CONFIDENCE') n8=np.array(n8) n8=n8.tolist() # n9 = f.get('LAND_COVERAGE') n9=np.array(n9) n9=n9.tolist() #Cria arquivo de saida fileName="/share/bramsrd/dist/BRAMS/data/FOCOS/{0:04d}{4:02d}{5:02d}/GOES-CAMS/f{0:04d}{1:03d}{2:02d}{3:02d}.v65.g16.filt".format(ano,djul,hora,min,mes,dia) #print("Opening %s" % fileName) f= open(fileName,"w+") #Escreve o cabeçalho de 5 linhas f.write("### Input "+fileIn+"\n") f.write("### \n") f.write("### Total de focos: {0:012d}\n".format(len(n3))) f.write("### Latitude, Longitude, ... \n") f.write("### deg, deg, ...\n") #Percorre os dados de entrada e escreve o arquivo de saída for i in range(len(n3)): #print("LAT: {0:6.2f} , LON: {1:6.2f}".format(n3[i],n2[i])) if n3[i]>-90 and n3[i]<-30.0: if n2[i]>-60.0 and n2[i]<15: #Determina o ponto em USGS para cada lat (n2) e lon (n3) for j in range(nLats): if n2[i]<=lat[j] and n2[i]>=lat[j+1]: pLat=j break for j in range(nLons): if n3[i]>=lon[j] and n3[i]<=lon[j+1]: pLon=j break #Pega o bioma desse ponto eco=imarray[pLat,pLon] #escreve os dados no arquivo de saída f.write("{0:9.4f} {1:9.4f} {2:9.4f} {3:9.4f} {4:9.4f} {5:9.4f} {8:6.3f} 0.0 {6:9.4f} 00.0 00 {7:03d} {9:03d}\n".format(n3[i],n2[i],n4[i],n5[i],n6[i],n7[i],n1[i],n8[i],fireArea[eco],eco)) npoints=npoints+1 #print("closing file %s " % fileName) f.close() print("Processados {0:6d} pontos para dados Goes-Cams!".format(npoints)) print("Escritos no arquivo "+fileName) print("") print("") if __name__ == '__main__': sys.exit(main(sys.argv))