diff --git a/README.md b/README.md index 882f33b..9c832bb 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,11 @@ In order to use the footprints code: -1) At the footprints project root (PATH-TO-PROJECT/FOOTPRINTS), set the environment variables: '. env.sh' +1) At the footprints project root (PATH-TO-PROJECT/FOOTPRINTS), +set the environment variables: '. env.sh' -2) Go back to the projects directory (/usr/local/soft) and set the PY3 python virtual environment: '. PY3/bin/activate' +2) Go back to the projects directory (/usr/local/soft) +and set the PY3 python virtual environment: '. PY3/bin/activate' -3) Finally, go to the footprints project and to the bin directory (/usr/local/soft/FOOTPRINTS/BIN) and run the code: 'python runFootprints.py (XML mission file)' or 'python runAll.py' (to run all missions described into the missions directory (/usr/local/soft/FOOTPRINTS/MISSIONS)) +3) Finally, go to the footprints project and to the +bin directory (PATH-TO-PROJECT/FOOTPRINTS/bin) and run the code: +'python runFootprints.py (XML_mission_file)' diff --git a/bin.py/runFootprints.py b/bin.py/runFootprints.py new file mode 100755 index 0000000..f7e8f3d --- /dev/null +++ b/bin.py/runFootprints.py @@ -0,0 +1,194 @@ +from lxml import etree +import pandas as pd +import time as myTime + +from tools import * + +def main(): + + # CHECK ARGS + checkArgs(2, 'Usage : python runFootprints.py (xml config file)') + + if os.path.isfile(os.getenv('FOOTPRINTSMISSIONS')+'/'+sys.argv[1]) == False: + print('[ERROR] '+sys.argv[1]+' does not exist') + sys.exit(2) + + start_time = myTime.time() + + # CONNECT TO DD SERVER + os.system('DDLogin biegun arnaud') + + # CREATE A DIRECTORY FOR THE NEW QUERY + MYDIR = os.getenv('FOOTPRINTSDATA')+'/'+(sys.argv[1]).split('.')[0]+'_'+datetime.now().strftime('%Y%m%dT%H%M%S') + os.system('mkdir '+MYDIR) + os.system('chmod ug+w ' + MYDIR) + + # CREATION OF A LOGGER + logFile = (sys.argv[1]).split('.')[0]+'_'+datetime.now().strftime('%Y%m%dT%H%M%S')+'.log' + logger = initLogger(MYDIR+'/'+logFile) + + # LOAD MISSION FILE + missionFile = os.getenv('FOOTPRINTSMISSIONS')+'/'+sys.argv[1] + satVI, satOrbParam, satOrbUnit, plasVI, indexVI, startTime, stopTime = loadMission(missionFile) + + # GET DD TIME BOUNDARIES + plasmaShift = 300 + plasDDStart, plasDDInter = getDDTimeBoundaries(startTime, stopTime, plasmaShift) + + satShift = 3600 + satDDStart, satDDInter = getDDTimeBoundaries(startTime, stopTime, satShift) + + indexShift = 300 + indexDDStart, indexDDInter = getDDTimeBoundaries(startTime, stopTime, indexShift) + + # PRINT MISSION INFOS + logger.info("Satellite : " + satVI) + logger.info("Plasma : " + plasVI) + logger.info("Sym index : " + indexVI) + logger.info("Start : " + startTime) + logger.info("Stop : " + stopTime) + logger.info("Plas DD Start : " + plasDDStart) + logger.info("Plas DD Inter : " + plasDDInter) + logger.info("Sat DD Start : " + satDDStart) + logger.info("Sat DD Inter : " + satDDInter) + + # GET DATA + logger.info('Loading data ...') + args = [plasVI,plasDDStart,plasDDInter,satVI,satOrbParam,satDDStart,satDDInter,indexVI,indexDDStart,indexDDInter] + getDataArgs = ' '.join(args) + os.system("GetData " + getDataArgs) + + if os.path.getsize('plasma.txt') == 0: + logger.error('Failed to load plasma data') + sys.exit(2) + elif os.path.getsize('sym.txt') == 0: + logger.error('Failed to load source orbit data') + sys.exit(2) + elif os.path.getsize('orbit.txt') == 0: + logger.error('Failed to load target orbit data') + sys.exit(2) + + os.system("mv plasma.txt "+MYDIR) + os.system("mv sym.txt "+MYDIR) + os.system("mv orbit.txt "+MYDIR) + + plasmaData = pd.read_csv(MYDIR+'/plasma.txt', dtype='S16,f4,f4,f4,f4,f4,f4,f4,f4') + plasmaData['Time'] = ddTime2Datetime(plasmaData['Time']) + + satData = pd.read_csv(MYDIR+'/orbit.txt', dtype='S16,f4,f4,f4') + satData['Time'] = ddTime2Datetime(satData['Time']) + + symData = pd.read_csv(MYDIR+'/sym.txt', dtype='S16,f4') + symData['Time'] = ddTime2Datetime(symData['Time']) + + # QUALITY FLAG FOR PLASMA DATA + qualFlag = np.isnan(plasmaData['vx']).astype(int) + qualData = {'Time': plasmaData['Time'],'QualityFlag': qualFlag} + qualFlagData = pd.DataFrame(data=qualData) + qualFlagData = qualFlagData.set_index('Time') + + # CLEAN PLASMA DATA + cleanedPlasmaData = plasmaData.interpolate().fillna(method='bfill') + + # BUILD THE MAIN DATAFRAME + myPlasma = cleanedPlasmaData.set_index('Time') + mySat = satData.set_index('Time') + mySym = symData.set_index('Time') + combined = myPlasma.join(mySat,how='outer').join(mySym,how='outer').join(qualFlagData,how='outer') + trunk = combined.truncate(before=startTime, after=stopTime) + df = trunk.groupby(trunk.index).first() + df = df.interpolate().fillna(method='backfill') + df = df.resample('5min') + df = df[['x', 'y', 'z', 'bx', 'by','bz','vx','vy','vz','ramP','symIndex','QualityFlag']] + df['symIndex'] = df['symIndex'].astype(int) + df['QualityFlag'] = df['QualityFlag'].astype(int) + + # NORMALIZE SPACE STEP TO EARTH RADIUS + if satOrbUnit != 're': + pc = loadPCFile(os.getenv('FOOTPRINTSCONFIG')+'/physical_constants.xml') + + if satOrbUnit == 'km': + df['x'] = df['x']/pc['RE'] + df['y'] = df['y']/pc['RE'] + df['z'] = df['z']/pc['RE'] + elif satOrbUnit == 'au': + df['x'] = df['x']*(pc['AU']/pc['RE']) + df['y'] = df['y']*(pc['AU']/pc['RE']) + df['z'] = df['z']*(pc['AU']/pc['RE']) + else: + logger.error('Unit for orbit data is unknown : '+satOrbUnit) + sys.exit(2) + + + # WRITE INPUTS FILE + prefix = satVI.split('_')[0] + '_' + startTime.split('-')[0] + inputsFile = MYDIR+'/'+prefix+'.txt' + writeInputFile(df,inputsFile) + + # COMPUTES FOOTPRINTS + fpFile = MYDIR+'/'+prefix+'_footprints.txt' + ncFile = MYDIR+'/'+prefix+'_footprints.nc' + logger.info('Computation of footprints ...') + os.system("Footprints " + inputsFile + " > " + fpFile) + os.system("NC " + fpFile + " " + ncFile) + + # END + stop_time = myTime.time() + interval = stop_time - start_time + logger.info('Program terminated in ' + str(interval) + ' seconds.') + + +############################## FUNCTIONS DEFINITIONS ################################## +def loadPCFile(pcFile): + myXml = etree.parse(pcFile) + au = ((myXml.xpath("/CONSTANTS/AU"))[0].text) + re = ((myXml.xpath("/CONSTANTS/RE"))[0].text) + + pc = { + 'AU': float(au), + 'RE': float(re), + } + + return pc + +def loadMission(missionFile): + myXml = etree.parse(missionFile) + satName = ((myXml.xpath("/MISSION/SATELLITE/NAME"))[0].text).lower() + satVI = ((myXml.xpath("/MISSION/SATELLITE/VI"))[0].text) + satOrbParam = ((myXml.xpath("/MISSION/SATELLITE/ORBIT"))[0].text) + satOrbUnit = ((myXml.xpath("/MISSION/SATELLITE/UNIT"))[0].text).lower() + plasVI = ((myXml.xpath("/MISSION/PLASMA/VI"))[0].text) + indexVI = ((myXml.xpath("/MISSION/INDEX/VI"))[0].text) + startTime = (myXml.xpath("/MISSION/START"))[0].text + stopTime = (myXml.xpath("/MISSION/STOP"))[0].text + + return satVI, satOrbParam, satOrbUnit, plasVI, indexVI, startTime, stopTime + +def writeInputFile(dataframe,inputsfilename): + with open(inputsfilename, 'w') as f: + + for i,j in enumerate(dataframe.index): + t = ' '.join([str(j.year),str(j.dayofyear),str(j.hour),str(j.minute)]) + x = '{:.2f}'.format(dataframe['x'][i]) + y = '{:.2f}'.format(dataframe['y'][i]) + z = '{:.2f}'.format(dataframe['z'][i]) + bx = '{:.2f}'.format(dataframe['bx'][i]) + by = '{:.2f}'.format(dataframe['by'][i]) + bz = '{:.2f}'.format(dataframe['bz'][i]) + vx = '{:.2f}'.format(dataframe['vx'][i]) + vy = '{:.2f}'.format(dataframe['vy'][i]) + vz = '{:.2f}'.format(dataframe['vz'][i]) + ramP = '{:.2f}'.format(dataframe['ramP'][i]) + sym = '{:.0f}'.format(dataframe['symIndex'][i]) + qF = '{:.0f}'.format(dataframe['QualityFlag'][i]) + params = [t,x,y,z,bx,by,bz,vx,vy,vz,ramP,sym,qF] + line = ' '.join(params)+'\n' + + f.write(line) + f.close() +############################### END OF FUNCTIONS DEFINITIONS ############################ + + + +if __name__ == '__main__': + main() diff --git a/bin.py/tools.py b/bin.py/tools.py new file mode 100755 index 0000000..78a42f1 --- /dev/null +++ b/bin.py/tools.py @@ -0,0 +1,177 @@ +######################## SYSTEM FUNCTIONS ############################## +import sys +import os + +def checkArgs(nbArgs, message): + if len(sys.argv) != nbArgs: + print(message) + sys.exit(2) +################## END OF SYSTEM FUNCTIONS ############################## + +######################## HANDLING ERRORS FUNCTIONS ###################### +import logging +from logging.handlers import RotatingFileHandler +# from logging.handlers import SMTPHandler + +def initLogger(logfile): + logger = logging.getLogger() + logger.setLevel(logging.INFO) + + fileFormatter = logging.Formatter('%(asctime)s :: [%(levelname)s] :: %(message)s') + fileHandler = RotatingFileHandler(logfile, 'a', 1000000, 1) + fileHandler.setLevel(logging.ERROR) + fileHandler.setFormatter(fileFormatter) + logger.addHandler(fileHandler) + + consoleFormatter = logging.Formatter('[%(levelname)s] %(message)s') + consoleHandler = logging.StreamHandler() + consoleHandler.setLevel(logging.INFO) + consoleHandler.setFormatter(consoleFormatter) + logger.addHandler(consoleHandler) + + # mailHandler = SMTPHandler( + # # Host et port + # 'webmail.irap.omp.eu', + # # From + # 'Arnaud.Biegun@irap.omp.eu', + # # To (liste) + # ['Arnaud.Biegun@irap.omp.eu'], + # # Sujet du message + # "Critical error when processing Solar Wind MHD model" + # ) + # mailHandler.setLevel(logging.DEBUG) + # logger.addHandler(mailHandler) + + return logger +######################## END OF HANDLING ERRORS ########################## + + +######################## TIME FUNCTIONS ############################## +from datetime import datetime, date, time, timedelta + +def time2ddtime(myDate): + fmt = "%Y-%m-%dT%H:%M:%S" + d = datetime.strptime(myDate, fmt) + dddoy = '{:03}'.format(int(d.strftime("%j"))-1) + hour = '{:02}'.format(d.hour) + minute = '{:02}'.format(d.minute) + second = '{:02}'.format(d.second) + + ddtime = str(d.year) + dddoy + hour + minute + second +'000' + + return ddtime + +def time2Unixtime(myDate): + fmt = "%Y-%m-%dT%H:%M:%S" + d = datetime.datetime.strptime(myDate, fmt) + + return d.strftime("%s") + +def time2J2000(myDate, j2000): + fmt = "%Y-%m-%dT%H:%M:%S" + d = int(datetime.datetime.strptime(myDate, fmt).strftime("%s")) + + return str(d-j2000) + + +def timeDelta(date1, date2): + fmt = "%Y-%m-%dT%H:%M:%S" + d1 = datetime.strptime(date1, fmt) + d2 = datetime.strptime(date2, fmt) + td = d2 - d1 + + intDays = td.days + intHours = td.seconds//3600 + intMinutes = (td.seconds//60)%60 + intSeconds = td.seconds-intHours*3600-intMinutes*60 + + return intDays, intHours, intMinutes, intSeconds + + +def DDTimeDelta(date1, date2): + interD, interH, interM, interS = timeDelta(date1, date2) + ddTimeDelta = '0000' + '{:03}'.format(interD) + '{:02}'.format(interH) + '{:02}'.format(interM) + '{:02}'.format(interS) + '000' + + return ddTimeDelta + +def shiftedDate(myDate, days=0, hours=0, minutes=0, seconds=0): + fmt = "%Y-%m-%dT%H:%M:%S" + d = datetime.strptime(myDate, fmt) + + myShiftedDate = d + timedelta(days=days, hours=hours, minutes=minutes, seconds=seconds) + + return myShiftedDate.strftime("%Y-%m-%dT%H:%M:%S") + +def getDDTimeBoundaries(startTime, stopTime, shift=0): + shiftedStart = shiftedDate(startTime, seconds=-shift) + shiftedStop = shiftedDate(stopTime, seconds=+shift) + DDStart = time2ddtime(shiftedStart) + DDInter = DDTimeDelta(shiftedStart, shiftedStop) + + return DDStart, DDInter + + +def ddTime2Datetime(ddTime): + myDatetime = [] + for t in ddTime: + year = int(t[0:4]) + day = int(t[4:7]) + 1 + myDate = datetime(year, 1, 1) + timedelta(day - 1) + + hour = int(t[7:9]) + minute = int(t[9:11]) + seconds = int(t[11:13]) + ms = int(t[13:16])*1000 + myTime = time(hour, minute,seconds,ms) + + myDatetime.append( datetime.combine(myDate, myTime) ) + + return myDatetime + + +def ddTime2Timestamp(ddTime): + newTime = [] + for t in ddTime: + year = t[0:4] + day = '{:03}'.format(int(t[4:7]) + 1) + other = t[7:] + newDD = year+day+other + newTime.append(newDD) + + myFormat = '%Y%j%H%M%S%f' + myTime = [ datetime.datetime.strptime( nti, myFormat ) for nti in newTime ] + timestamps = [ time.mktime(mti.timetuple()) for mti in myTime ] + + return timestamps + +def getDDTimeFromJ2000(j2000Time, j2000): + + fmt = "%Y-%m-%d %H:%M:%S" + fromJ2000 = [ datetime.datetime.strptime(str(datetime.datetime.fromtimestamp(ti+j2000)), fmt) for ti in j2000Time ] + dd = [] + for ti in fromJ2000: + dddoy = '{:03}'.format(int(ti.strftime("%j"))-1) + hour = '{:02}'.format(ti.hour) + minute = '{:02}'.format(ti.minute) + second = '{:02}'.format(ti.second) + + ddtime = str(ti.year) + dddoy + hour + minute + second +'000' + + dd.append(ddtime) + + return dd + +################## END OF TIME FUNCTIONS ############################## + + +######################## DATA FUNCTIONS ############################## +import numpy as np + +def cleanData(data): + + # Clean NAN values + mask = np.isnan(data) + data[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask]) + + return data +################## END OF DATA FUNCTIONS ############################## \ No newline at end of file diff --git a/config/DD.res b/config/DD.res new file mode 100644 index 0000000..e8a9739 --- /dev/null +++ b/config/DD.res @@ -0,0 +1,2 @@ +manunja.irap.omp.eu +5600 \ No newline at end of file diff --git a/config/physical_constants.xml b/config/physical_constants.xml new file mode 100644 index 0000000..5714529 --- /dev/null +++ b/config/physical_constants.xml @@ -0,0 +1,5 @@ + + + 149597871 + 6371 + \ No newline at end of file diff --git a/env.sh b/env.sh index df82260..a44c61d 100644 --- a/env.sh +++ b/env.sh @@ -11,11 +11,11 @@ export TINYDIR=/usr/local/tinyxml export BOOSTDIR=/usr/local/boost #------- DO NOT EDIT BELOW -------- -export FOOTPRINTSMISSIONS=$FOOTPRINTSROOT/MISSIONS -export FOOTPRINTSCONFIG=$FOOTPRINTSROOT/CONFIG +export FOOTPRINTSMISSIONS=${FOOTPRINTSROOT}/missions +export FOOTPRINTSCONFIG=$FOOTPRINTSROOT/config export LD_LIBRARY_PATH=$GEOPACKDIR/lib:$BOOSTDIR/stage/lib:$DDDIR/lib:$TINYDIR/lib -export PATH=/bin:/usr/bin:/usr/local/bin:./:$FOOTPRINTSROOT/BIN:$DDDIR/bin - +export PATH=/bin:/usr/bin:/usr/local/bin:./:$FOOTPRINTSROOT/bin:$DDDIR/bin +export DDPATH=${FOOTPRINTSCONFIG} # DATA ROOT -export FOOTPRINTSDATA=$FOOTPRINTSROOT/DATA -export FOOTPRINTSFILES=$FOOTPRINTSDATA/FOOTPRINTS \ No newline at end of file +export FOOTPRINTSDATA=$FOOTPRINTSROOT/data +export FOOTPRINTSFILES=$FOOTPRINTSDATA/footprints \ No newline at end of file diff --git a/missions/C1_2016.xml b/missions/C1_2016.xml new file mode 100644 index 0000000..6842eb4 --- /dev/null +++ b/missions/C1_2016.xml @@ -0,0 +1,25 @@ + + + + CLUSTER_1 + clust1_orb_all + Orb + 60 + Re + + + OMNI + omni_5min_all + V + RamP + Vel + B + + + SYM_H + ground_based_asy + SYM + + 2016-01-01T00:00:00 + 2016-03-15T23:00:00 + diff --git a/src/Makefile b/src/Makefile index 0827680..95de3c0 100644 --- a/src/Makefile +++ b/src/Makefile @@ -16,10 +16,9 @@ EXE = GetData Footprints NC # ******************************** LINKING ********************************* all: $(EXE) install -d -m a+rx,ug+w ../bin - # mv Inputs ../bin - mv Footprints ../bin - mv NC ../bin - + mv ${EXE} ../bin + install -d -m a+rx,ug+w ../data + NC: makeNcFile.o ctools.o gcc -o $@ $^ -lnetcdf -lm -- libgit2 0.21.2