CDPP / FOOTPRINTS

Commit 68dc3ad5bd1f546c4a966d8a76b49558f85d612f

Authored by Elena.Budnik
1 parent 7c02e6a3
Exists in master

all project

Showing 8 changed files with 419 additions and 13 deletions   Show diff stats
Inline Side-by-side
README.md
  Diff comments   View file @68dc3ad
1 1 In order to use the footprints code:
2 2  
3   -1) At the footprints project root (PATH-TO-PROJECT/FOOTPRINTS), set the environment variables: '. env.sh'
  3 +1) At the footprints project root (PATH-TO-PROJECT/FOOTPRINTS),
  4 +set the environment variables: '. env.sh'
4 5  
5   -2) Go back to the projects directory (/usr/local/soft) and set the PY3 python virtual environment: '. PY3/bin/activate'
  6 +2) Go back to the projects directory (/usr/local/soft)
  7 +and set the PY3 python virtual environment: '. PY3/bin/activate'
6 8  
7   -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))
  9 +3) Finally, go to the footprints project and to the
  10 +bin directory (PATH-TO-PROJECT/FOOTPRINTS/bin) and run the code:
  11 +'python runFootprints.py (XML_mission_file)'
... ...
bin.py/runFootprints.py 0 → 100755
  Diff comments   View file @68dc3ad
... ... @@ -0,0 +1,194 @@
  1 +from lxml import etree
  2 +import pandas as pd
  3 +import time as myTime
  4 +
  5 +from tools import *
  6 +
  7 +def main():
  8 +
  9 + # CHECK ARGS
  10 + checkArgs(2, 'Usage : python runFootprints.py (xml config file)')
  11 +
  12 + if os.path.isfile(os.getenv('FOOTPRINTSMISSIONS')+'/'+sys.argv[1]) == False:
  13 + print('[ERROR] '+sys.argv[1]+' does not exist')
  14 + sys.exit(2)
  15 +
  16 + start_time = myTime.time()
  17 +
  18 + # CONNECT TO DD SERVER
  19 + os.system('DDLogin biegun arnaud')
  20 +
  21 + # CREATE A DIRECTORY FOR THE NEW QUERY
  22 + MYDIR = os.getenv('FOOTPRINTSDATA')+'/'+(sys.argv[1]).split('.')[0]+'_'+datetime.now().strftime('%Y%m%dT%H%M%S')
  23 + os.system('mkdir '+MYDIR)
  24 + os.system('chmod ug+w ' + MYDIR)
  25 +
  26 + # CREATION OF A LOGGER
  27 + logFile = (sys.argv[1]).split('.')[0]+'_'+datetime.now().strftime('%Y%m%dT%H%M%S')+'.log'
  28 + logger = initLogger(MYDIR+'/'+logFile)
  29 +
  30 + # LOAD MISSION FILE
  31 + missionFile = os.getenv('FOOTPRINTSMISSIONS')+'/'+sys.argv[1]
  32 + satVI, satOrbParam, satOrbUnit, plasVI, indexVI, startTime, stopTime = loadMission(missionFile)
  33 +
  34 + # GET DD TIME BOUNDARIES
  35 + plasmaShift = 300
  36 + plasDDStart, plasDDInter = getDDTimeBoundaries(startTime, stopTime, plasmaShift)
  37 +
  38 + satShift = 3600
  39 + satDDStart, satDDInter = getDDTimeBoundaries(startTime, stopTime, satShift)
  40 +
  41 + indexShift = 300
  42 + indexDDStart, indexDDInter = getDDTimeBoundaries(startTime, stopTime, indexShift)
  43 +
  44 + # PRINT MISSION INFOS
  45 + logger.info("Satellite : " + satVI)
  46 + logger.info("Plasma : " + plasVI)
  47 + logger.info("Sym index : " + indexVI)
  48 + logger.info("Start : " + startTime)
  49 + logger.info("Stop : " + stopTime)
  50 + logger.info("Plas DD Start : " + plasDDStart)
  51 + logger.info("Plas DD Inter : " + plasDDInter)
  52 + logger.info("Sat DD Start : " + satDDStart)
  53 + logger.info("Sat DD Inter : " + satDDInter)
  54 +
  55 + # GET DATA
  56 + logger.info('Loading data ...')
  57 + args = [plasVI,plasDDStart,plasDDInter,satVI,satOrbParam,satDDStart,satDDInter,indexVI,indexDDStart,indexDDInter]
  58 + getDataArgs = ' '.join(args)
  59 + os.system("GetData " + getDataArgs)
  60 +
  61 + if os.path.getsize('plasma.txt') == 0:
  62 + logger.error('Failed to load plasma data')
  63 + sys.exit(2)
  64 + elif os.path.getsize('sym.txt') == 0:
  65 + logger.error('Failed to load source orbit data')
  66 + sys.exit(2)
  67 + elif os.path.getsize('orbit.txt') == 0:
  68 + logger.error('Failed to load target orbit data')
  69 + sys.exit(2)
  70 +
  71 + os.system("mv plasma.txt "+MYDIR)
  72 + os.system("mv sym.txt "+MYDIR)
  73 + os.system("mv orbit.txt "+MYDIR)
  74 +
  75 + plasmaData = pd.read_csv(MYDIR+'/plasma.txt', dtype='S16,f4,f4,f4,f4,f4,f4,f4,f4')
  76 + plasmaData['Time'] = ddTime2Datetime(plasmaData['Time'])
  77 +
  78 + satData = pd.read_csv(MYDIR+'/orbit.txt', dtype='S16,f4,f4,f4')
  79 + satData['Time'] = ddTime2Datetime(satData['Time'])
  80 +
  81 + symData = pd.read_csv(MYDIR+'/sym.txt', dtype='S16,f4')
  82 + symData['Time'] = ddTime2Datetime(symData['Time'])
  83 +
  84 + # QUALITY FLAG FOR PLASMA DATA
  85 + qualFlag = np.isnan(plasmaData['vx']).astype(int)
  86 + qualData = {'Time': plasmaData['Time'],'QualityFlag': qualFlag}
  87 + qualFlagData = pd.DataFrame(data=qualData)
  88 + qualFlagData = qualFlagData.set_index('Time')
  89 +
  90 + # CLEAN PLASMA DATA
  91 + cleanedPlasmaData = plasmaData.interpolate().fillna(method='bfill')
  92 +
  93 + # BUILD THE MAIN DATAFRAME
  94 + myPlasma = cleanedPlasmaData.set_index('Time')
  95 + mySat = satData.set_index('Time')
  96 + mySym = symData.set_index('Time')
  97 + combined = myPlasma.join(mySat,how='outer').join(mySym,how='outer').join(qualFlagData,how='outer')
  98 + trunk = combined.truncate(before=startTime, after=stopTime)
  99 + df = trunk.groupby(trunk.index).first()
  100 + df = df.interpolate().fillna(method='backfill')
  101 + df = df.resample('5min')
  102 + df = df[['x', 'y', 'z', 'bx', 'by','bz','vx','vy','vz','ramP','symIndex','QualityFlag']]
  103 + df['symIndex'] = df['symIndex'].astype(int)
  104 + df['QualityFlag'] = df['QualityFlag'].astype(int)
  105 +
  106 + # NORMALIZE SPACE STEP TO EARTH RADIUS
  107 + if satOrbUnit != 're':
  108 + pc = loadPCFile(os.getenv('FOOTPRINTSCONFIG')+'/physical_constants.xml')
  109 +
  110 + if satOrbUnit == 'km':
  111 + df['x'] = df['x']/pc['RE']
  112 + df['y'] = df['y']/pc['RE']
  113 + df['z'] = df['z']/pc['RE']
  114 + elif satOrbUnit == 'au':
  115 + df['x'] = df['x']*(pc['AU']/pc['RE'])
  116 + df['y'] = df['y']*(pc['AU']/pc['RE'])
  117 + df['z'] = df['z']*(pc['AU']/pc['RE'])
  118 + else:
  119 + logger.error('Unit for orbit data is unknown : '+satOrbUnit)
  120 + sys.exit(2)
  121 +
  122 +
  123 + # WRITE INPUTS FILE
  124 + prefix = satVI.split('_')[0] + '_' + startTime.split('-')[0]
  125 + inputsFile = MYDIR+'/'+prefix+'.txt'
  126 + writeInputFile(df,inputsFile)
  127 +
  128 + # COMPUTES FOOTPRINTS
  129 + fpFile = MYDIR+'/'+prefix+'_footprints.txt'
  130 + ncFile = MYDIR+'/'+prefix+'_footprints.nc'
  131 + logger.info('Computation of footprints ...')
  132 + os.system("Footprints " + inputsFile + " > " + fpFile)
  133 + os.system("NC " + fpFile + " " + ncFile)
  134 +
  135 + # END
  136 + stop_time = myTime.time()
  137 + interval = stop_time - start_time
  138 + logger.info('Program terminated in ' + str(interval) + ' seconds.')
  139 +
  140 +
  141 +############################## FUNCTIONS DEFINITIONS ##################################
  142 +def loadPCFile(pcFile):
  143 + myXml = etree.parse(pcFile)
  144 + au = ((myXml.xpath("/CONSTANTS/AU"))[0].text)
  145 + re = ((myXml.xpath("/CONSTANTS/RE"))[0].text)
  146 +
  147 + pc = {
  148 + 'AU': float(au),
  149 + 'RE': float(re),
  150 + }
  151 +
  152 + return pc
  153 +
  154 +def loadMission(missionFile):
  155 + myXml = etree.parse(missionFile)
  156 + satName = ((myXml.xpath("/MISSION/SATELLITE/NAME"))[0].text).lower()
  157 + satVI = ((myXml.xpath("/MISSION/SATELLITE/VI"))[0].text)
  158 + satOrbParam = ((myXml.xpath("/MISSION/SATELLITE/ORBIT"))[0].text)
  159 + satOrbUnit = ((myXml.xpath("/MISSION/SATELLITE/UNIT"))[0].text).lower()
  160 + plasVI = ((myXml.xpath("/MISSION/PLASMA/VI"))[0].text)
  161 + indexVI = ((myXml.xpath("/MISSION/INDEX/VI"))[0].text)
  162 + startTime = (myXml.xpath("/MISSION/START"))[0].text
  163 + stopTime = (myXml.xpath("/MISSION/STOP"))[0].text
  164 +
  165 + return satVI, satOrbParam, satOrbUnit, plasVI, indexVI, startTime, stopTime
  166 +
  167 +def writeInputFile(dataframe,inputsfilename):
  168 + with open(inputsfilename, 'w') as f:
  169 +
  170 + for i,j in enumerate(dataframe.index):
  171 + t = ' '.join([str(j.year),str(j.dayofyear),str(j.hour),str(j.minute)])
  172 + x = '{:.2f}'.format(dataframe['x'][i])
  173 + y = '{:.2f}'.format(dataframe['y'][i])
  174 + z = '{:.2f}'.format(dataframe['z'][i])
  175 + bx = '{:.2f}'.format(dataframe['bx'][i])
  176 + by = '{:.2f}'.format(dataframe['by'][i])
  177 + bz = '{:.2f}'.format(dataframe['bz'][i])
  178 + vx = '{:.2f}'.format(dataframe['vx'][i])
  179 + vy = '{:.2f}'.format(dataframe['vy'][i])
  180 + vz = '{:.2f}'.format(dataframe['vz'][i])
  181 + ramP = '{:.2f}'.format(dataframe['ramP'][i])
  182 + sym = '{:.0f}'.format(dataframe['symIndex'][i])
  183 + qF = '{:.0f}'.format(dataframe['QualityFlag'][i])
  184 + params = [t,x,y,z,bx,by,bz,vx,vy,vz,ramP,sym,qF]
  185 + line = ' '.join(params)+'\n'
  186 +
  187 + f.write(line)
  188 + f.close()
  189 +############################### END OF FUNCTIONS DEFINITIONS ############################
  190 +
  191 +
  192 +
  193 +if __name__ == '__main__':
  194 + main()
... ...
bin.py/tools.py 0 → 100755
  Diff comments   View file @68dc3ad
... ... @@ -0,0 +1,177 @@
  1 +######################## SYSTEM FUNCTIONS ##############################
  2 +import sys
  3 +import os
  4 +
  5 +def checkArgs(nbArgs, message):
  6 + if len(sys.argv) != nbArgs:
  7 + print(message)
  8 + sys.exit(2)
  9 +################## END OF SYSTEM FUNCTIONS ##############################
  10 +
  11 +######################## HANDLING ERRORS FUNCTIONS ######################
  12 +import logging
  13 +from logging.handlers import RotatingFileHandler
  14 +# from logging.handlers import SMTPHandler
  15 +
  16 +def initLogger(logfile):
  17 + logger = logging.getLogger()
  18 + logger.setLevel(logging.INFO)
  19 +
  20 + fileFormatter = logging.Formatter('%(asctime)s :: [%(levelname)s] :: %(message)s')
  21 + fileHandler = RotatingFileHandler(logfile, 'a', 1000000, 1)
  22 + fileHandler.setLevel(logging.ERROR)
  23 + fileHandler.setFormatter(fileFormatter)
  24 + logger.addHandler(fileHandler)
  25 +
  26 + consoleFormatter = logging.Formatter('[%(levelname)s] %(message)s')
  27 + consoleHandler = logging.StreamHandler()
  28 + consoleHandler.setLevel(logging.INFO)
  29 + consoleHandler.setFormatter(consoleFormatter)
  30 + logger.addHandler(consoleHandler)
  31 +
  32 + # mailHandler = SMTPHandler(
  33 + # # Host et port
  34 + # 'webmail.irap.omp.eu',
  35 + # # From
  36 + # 'Arnaud.Biegun@irap.omp.eu',
  37 + # # To (liste)
  38 + # ['Arnaud.Biegun@irap.omp.eu'],
  39 + # # Sujet du message
  40 + # "Critical error when processing Solar Wind MHD model"
  41 + # )
  42 + # mailHandler.setLevel(logging.DEBUG)
  43 + # logger.addHandler(mailHandler)
  44 +
  45 + return logger
  46 +######################## END OF HANDLING ERRORS ##########################
  47 +
  48 +
  49 +######################## TIME FUNCTIONS ##############################
  50 +from datetime import datetime, date, time, timedelta
  51 +
  52 +def time2ddtime(myDate):
  53 + fmt = "%Y-%m-%dT%H:%M:%S"
  54 + d = datetime.strptime(myDate, fmt)
  55 + dddoy = '{:03}'.format(int(d.strftime("%j"))-1)
  56 + hour = '{:02}'.format(d.hour)
  57 + minute = '{:02}'.format(d.minute)
  58 + second = '{:02}'.format(d.second)
  59 +
  60 + ddtime = str(d.year) + dddoy + hour + minute + second +'000'
  61 +
  62 + return ddtime
  63 +
  64 +def time2Unixtime(myDate):
  65 + fmt = "%Y-%m-%dT%H:%M:%S"
  66 + d = datetime.datetime.strptime(myDate, fmt)
  67 +
  68 + return d.strftime("%s")
  69 +
  70 +def time2J2000(myDate, j2000):
  71 + fmt = "%Y-%m-%dT%H:%M:%S"
  72 + d = int(datetime.datetime.strptime(myDate, fmt).strftime("%s"))
  73 +
  74 + return str(d-j2000)
  75 +
  76 +
  77 +def timeDelta(date1, date2):
  78 + fmt = "%Y-%m-%dT%H:%M:%S"
  79 + d1 = datetime.strptime(date1, fmt)
  80 + d2 = datetime.strptime(date2, fmt)
  81 + td = d2 - d1
  82 +
  83 + intDays = td.days
  84 + intHours = td.seconds//3600
  85 + intMinutes = (td.seconds//60)%60
  86 + intSeconds = td.seconds-intHours*3600-intMinutes*60
  87 +
  88 + return intDays, intHours, intMinutes, intSeconds
  89 +
  90 +
  91 +def DDTimeDelta(date1, date2):
  92 + interD, interH, interM, interS = timeDelta(date1, date2)
  93 + ddTimeDelta = '0000' + '{:03}'.format(interD) + '{:02}'.format(interH) + '{:02}'.format(interM) + '{:02}'.format(interS) + '000'
  94 +
  95 + return ddTimeDelta
  96 +
  97 +def shiftedDate(myDate, days=0, hours=0, minutes=0, seconds=0):
  98 + fmt = "%Y-%m-%dT%H:%M:%S"
  99 + d = datetime.strptime(myDate, fmt)
  100 +
  101 + myShiftedDate = d + timedelta(days=days, hours=hours, minutes=minutes, seconds=seconds)
  102 +
  103 + return myShiftedDate.strftime("%Y-%m-%dT%H:%M:%S")
  104 +
  105 +def getDDTimeBoundaries(startTime, stopTime, shift=0):
  106 + shiftedStart = shiftedDate(startTime, seconds=-shift)
  107 + shiftedStop = shiftedDate(stopTime, seconds=+shift)
  108 + DDStart = time2ddtime(shiftedStart)
  109 + DDInter = DDTimeDelta(shiftedStart, shiftedStop)
  110 +
  111 + return DDStart, DDInter
  112 +
  113 +
  114 +def ddTime2Datetime(ddTime):
  115 + myDatetime = []
  116 + for t in ddTime:
  117 + year = int(t[0:4])
  118 + day = int(t[4:7]) + 1
  119 + myDate = datetime(year, 1, 1) + timedelta(day - 1)
  120 +
  121 + hour = int(t[7:9])
  122 + minute = int(t[9:11])
  123 + seconds = int(t[11:13])
  124 + ms = int(t[13:16])*1000
  125 + myTime = time(hour, minute,seconds,ms)
  126 +
  127 + myDatetime.append( datetime.combine(myDate, myTime) )
  128 +
  129 + return myDatetime
  130 +
  131 +
  132 +def ddTime2Timestamp(ddTime):
  133 + newTime = []
  134 + for t in ddTime:
  135 + year = t[0:4]
  136 + day = '{:03}'.format(int(t[4:7]) + 1)
  137 + other = t[7:]
  138 + newDD = year+day+other
  139 + newTime.append(newDD)
  140 +
  141 + myFormat = '%Y%j%H%M%S%f'
  142 + myTime = [ datetime.datetime.strptime( nti, myFormat ) for nti in newTime ]
  143 + timestamps = [ time.mktime(mti.timetuple()) for mti in myTime ]
  144 +
  145 + return timestamps
  146 +
  147 +def getDDTimeFromJ2000(j2000Time, j2000):
  148 +
  149 + fmt = "%Y-%m-%d %H:%M:%S"
  150 + fromJ2000 = [ datetime.datetime.strptime(str(datetime.datetime.fromtimestamp(ti+j2000)), fmt) for ti in j2000Time ]
  151 + dd = []
  152 + for ti in fromJ2000:
  153 + dddoy = '{:03}'.format(int(ti.strftime("%j"))-1)
  154 + hour = '{:02}'.format(ti.hour)
  155 + minute = '{:02}'.format(ti.minute)
  156 + second = '{:02}'.format(ti.second)
  157 +
  158 + ddtime = str(ti.year) + dddoy + hour + minute + second +'000'
  159 +
  160 + dd.append(ddtime)
  161 +
  162 + return dd
  163 +
  164 +################## END OF TIME FUNCTIONS ##############################
  165 +
  166 +
  167 +######################## DATA FUNCTIONS ##############################
  168 +import numpy as np
  169 +
  170 +def cleanData(data):
  171 +
  172 + # Clean NAN values
  173 + mask = np.isnan(data)
  174 + data[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), data[~mask])
  175 +
  176 + return data
  177 +################## END OF DATA FUNCTIONS ##############################
0 178 \ No newline at end of file
... ...
config/DD.res 0 → 100644
  Diff comments   View file @68dc3ad
... ... @@ -0,0 +1,2 @@
  1 +manunja.irap.omp.eu
  2 +5600
0 3 \ No newline at end of file
... ...
config/physical_constants.xml 0 → 100644
  Diff comments   View file @68dc3ad
... ... @@ -0,0 +1,5 @@
  1 +<?xml version="1.0" ?>
  2 +<CONSTANTS>
  3 + <AU>149597871</AU>
  4 + <RE>6371</RE>
  5 +</CONSTANTS>
0 6 \ No newline at end of file
... ...
env.sh
  Diff comments   View file @68dc3ad
... ... @@ -11,11 +11,11 @@ export TINYDIR=/usr/local/tinyxml
11 11 export BOOSTDIR=/usr/local/boost
12 12  
13 13 #------- DO NOT EDIT BELOW --------
14   -export FOOTPRINTSMISSIONS=$FOOTPRINTSROOT/MISSIONS
15   -export FOOTPRINTSCONFIG=$FOOTPRINTSROOT/CONFIG
  14 +export FOOTPRINTSMISSIONS=${FOOTPRINTSROOT}/missions
  15 +export FOOTPRINTSCONFIG=$FOOTPRINTSROOT/config
16 16 export LD_LIBRARY_PATH=$GEOPACKDIR/lib:$BOOSTDIR/stage/lib:$DDDIR/lib:$TINYDIR/lib
17   -export PATH=/bin:/usr/bin:/usr/local/bin:./:$FOOTPRINTSROOT/BIN:$DDDIR/bin
18   -
  17 +export PATH=/bin:/usr/bin:/usr/local/bin:./:$FOOTPRINTSROOT/bin:$DDDIR/bin
  18 +export DDPATH=${FOOTPRINTSCONFIG}
19 19 # DATA ROOT
20   -export FOOTPRINTSDATA=$FOOTPRINTSROOT/DATA
21   -export FOOTPRINTSFILES=$FOOTPRINTSDATA/FOOTPRINTS
22 20 \ No newline at end of file
  21 +export FOOTPRINTSDATA=$FOOTPRINTSROOT/data
  22 +export FOOTPRINTSFILES=$FOOTPRINTSDATA/footprints
23 23 \ No newline at end of file
... ...
missions/C1_2016.xml 0 → 100644
  Diff comments   View file @68dc3ad
... ... @@ -0,0 +1,25 @@
  1 +<?xml version="1.0" ?>
  2 +<MISSION>
  3 + <SATELLITE>
  4 + <NAME>CLUSTER_1</NAME>
  5 + <VI>clust1_orb_all</VI>
  6 + <ORBIT>Orb</ORBIT>
  7 + <RES>60</RES>
  8 + <UNIT>Re</UNIT>
  9 + </SATELLITE>
  10 + <PLASMA>
  11 + <NAME>OMNI</NAME>
  12 + <VI>omni_5min_all</VI>
  13 + <VNORM>V</VNORM>
  14 + <RAMP>RamP</RAMP>
  15 + <VEL>Vel</VEL>
  16 + <MAG>B</MAG>
  17 + </PLASMA>
  18 + <INDEX>
  19 + <NAME>SYM_H</NAME>
  20 + <VI>ground_based_asy</VI>
  21 + <SYM>SYM</SYM>
  22 + </INDEX>
  23 + <START>2016-01-01T00:00:00</START>
  24 + <STOP>2016-03-15T23:00:00</STOP>
  25 +</MISSION>
... ...
src/Makefile
  Diff comments   View file @68dc3ad
... ... @@ -16,10 +16,9 @@ EXE = GetData Footprints NC
16 16 # ******************************** LINKING *********************************
17 17 all: $(EXE)
18 18 install -d -m a+rx,ug+w ../bin
19   - # mv Inputs ../bin
20   - mv Footprints ../bin
21   - mv NC ../bin
22   -
  19 + mv ${EXE} ../bin
  20 + install -d -m a+rx,ug+w ../data
  21 +
23 22 NC: makeNcFile.o ctools.o
24 23 gcc -o $@ $^ -lnetcdf -lm
25 24  
... ...