155 lines
5.8 KiB
Python
155 lines
5.8 KiB
Python
#!/usr/bin/env python2
|
||
# -*- coding: UTF-8 -*-
|
||
#
|
||
# INPUTS:
|
||
# OUTPUTS：download files ****.py
|
||
import numpy as np
|
||
import scipy.io
|
||
import datetime
|
||
import os
|
||
import re
|
||
import struct
|
||
|
||
def aps_era_ECMWF_download(era_workdir,lonmin,lonmax,latmin,latmax,UTC,SLC1,SLC2):
|
||
|
||
UTC_sat = UTC
|
||
SLC_list=[int(SLC1),int(SLC2)]
|
||
#%% Compute based on Satellite pass which ERA-I times will be used
|
||
#find two closest times with respect the the 6 hr ERA-I data
|
||
timelist_ERA = ['0000','0600' ,'1200','1800' , '0000']
|
||
time = int(UTC_sat[0:2]) + float(UTC_sat[3:5])/60
|
||
#print(time)
|
||
t_before = int(np.floor(time/6))
|
||
t_after = int(np.ceil(time/6))
|
||
#the faction it is closer towards the other date.
|
||
f_after = (time - 6*t_before)/(6*t_after - 6*t_before)
|
||
f_before = 1-f_after
|
||
|
||
# the time stamp of the closest two ERA acquisitions
|
||
time1 = str(timelist_ERA[t_before])
|
||
time2 = str(timelist_ERA[t_after])
|
||
|
||
#download era
|
||
S = str(latmin - 2)
|
||
N = str(latmax + 2)
|
||
W = str(lonmin - 2)
|
||
E = str(lonmax + 2)
|
||
wheatherregstr = N+'/'+W+'/'+S+'/'+E; # N/W/S/E
|
||
|
||
for date in SLC_list:
|
||
date2=date
|
||
if t_after==4:
|
||
date=str(date)
|
||
date2=((datetime.datetime.strptime(date,'%Y%m%d')+datetime.timedelta(days=1)).date()).strftime('%Y%m%d')
|
||
|
||
os.chdir(era_workdir)
|
||
date_path=os.path.join(era_workdir,str(date))
|
||
if not(os.path.exists(date_path)):
|
||
os.makedirs(date_path)
|
||
os.chdir(date_path)
|
||
|
||
for kk in range(0,2):
|
||
if kk==0:
|
||
file = 'ggap'+str(date)+time1+'.nc'
|
||
date=str(date)
|
||
year = date[0:4]
|
||
month = date[4:6]
|
||
day = date[6:8]
|
||
datestr = year+'-'+month+'-'+day
|
||
htime=time1[0:2]
|
||
print(datestr,htime)
|
||
timing=str(date)+time1
|
||
pythonsc_path = timing+'.py'
|
||
with open(pythonsc_path,"wb") as g:
|
||
|
||
g.write('#!/usr/bin/env python\n');
|
||
g.write('\n');
|
||
g.write('from ecmwfapi import ECMWFDataServer\n');
|
||
g.write('\n');
|
||
g.write('# To run this example, you need an API key\n');
|
||
g.write('# available from https://api.ecmwf.int/v1/key/\n');
|
||
g.write('\n');
|
||
g.write('server = ECMWFDataServer()\n');
|
||
g.write('\n');
|
||
g.write('server.retrieve({\n');
|
||
g.write(' \'dataset\' : "interim",\n');
|
||
g.write(' \'class\' : "ei",\n');
|
||
g.write(' \'type\' : "an",\n');
|
||
g.write(' \'stream\' : "oper",\n');
|
||
g.write(' \'levtype\' : "pl",\n');
|
||
g.write(' \'levelist\': "1/2/3/5/7/10/20/30/50/70/100/125/150/175/200/225/250/300/350/400/450/500/550/600/650/700/750/775/800/825/850/875/900/925/950/975/1000",\n');
|
||
g.write(' \'param\' : "129.128/130.128/157.128/246.128",\n');
|
||
g.write(' \'date\' : "%s",\n'%(datestr));
|
||
g.write(' \'time\' : "%s",\n'%(htime));
|
||
g.write(' \'step\' : "0",\n');
|
||
g.write(' \'format\' : "netcdf",\n');
|
||
g.write(' \'resol\' : "auto",\n');
|
||
g.write(' \'area\' : "%s",\n'%(wheatherregstr));
|
||
g.write(' \'grid\' : "0.125/0.125",\n');
|
||
g.write(' \'target\' : "%s",\n'%(file ));
|
||
g.write(' })\n');
|
||
|
||
if kk==1:
|
||
file = 'ggap'+str(date2)+time2+'.nc'
|
||
date2=str(date2)
|
||
year = date2[0:4]
|
||
month = date2[4:6]
|
||
day = date2[6:8]
|
||
datestr = year+'-'+month+'-'+day
|
||
htime=time2[0:2]
|
||
print(datestr,htime)
|
||
timing=str(date2)+time2
|
||
pythonsc_path = timing+'.py'
|
||
with open(pythonsc_path,"wb") as g:
|
||
|
||
g.write('#!/usr/bin/env python\n');
|
||
g.write('\n');
|
||
g.write('from ecmwfapi import ECMWFDataServer\n');
|
||
g.write('\n');
|
||
g.write('# To run this example, you need an API key\n');
|
||
g.write('# available from https://api.ecmwf.int/v1/key/\n');
|
||
g.write('\n');
|
||
g.write('server = ECMWFDataServer()\n');
|
||
g.write('\n');
|
||
g.write('server.retrieve({\n');
|
||
g.write(' \'dataset\' : "interim",\n');
|
||
g.write(' \'class\' : "ei",\n');
|
||
g.write(' \'type\' : "an",\n');
|
||
g.write(' \'stream\' : "oper",\n');
|
||
g.write(' \'levtype\' : "pl",\n');
|
||
g.write(' \'levelist\': "1/2/3/5/7/10/20/30/50/70/100/125/150/175/200/225/250/300/350/400/450/500/550/600/650/700/750/775/800/825/850/875/900/925/950/975/1000",\n');
|
||
g.write(' \'param\' : "129.128/130.128/157.128/246.128",\n');
|
||
g.write(' \'date\' : "%s",\n'%(datestr));
|
||
g.write(' \'time\' : "%s",\n'%(htime));
|
||
g.write(' \'step\' : "0",\n');
|
||
g.write(' \'format\' : "netcdf",\n');
|
||
g.write(' \'resol\' : "auto",\n');
|
||
g.write(' \'area\' : "%s",\n'%(wheatherregstr));
|
||
g.write(' \'grid\' : "0.125/0.125",\n');
|
||
g.write(' \'target\' : "%s",\n'%(file ));
|
||
g.write(' })\n');
|
||
|
||
if __name__=="__main__":
|
||
import sys
|
||
if (len(sys.argv)<8 or len(sys.argv)>9):
|
||
print('[1]era_workdir (input) The path of era files')
|
||
print('[2]lonmin (input) The min longitude. Make this slightly bigger than your InSAR region ')
|
||
print('[3]lonmax (input) The max longitude. Make this slightly bigger than your InSAR region' )
|
||
print('[4]latmin (input) The min latitude. Make this slightly bigger than your InSAR region')
|
||
print('[5]latmax (input) The max latitude. Make this slightly bigger than your InSAR region')
|
||
print('[6]UTC (input) The UTC time of radar image :14:44')
|
||
print('[7]SLC1 (input) The radar image :20160829')
|
||
print('[8]SLC2 (input) The radar image :20161001')
|
||
|
||
sys.exit()
|
||
era_workdir=sys.argv[1]
|
||
lonmin=float(sys.argv[2])
|
||
lonmax=float(sys.argv[3])
|
||
latmin=float(sys.argv[4])
|
||
latmax=float(sys.argv[5])
|
||
UTC=sys.argv[6]
|
||
SLC1=sys.argv[7]
|
||
SLC2=sys.argv[8]
|
||
|
||
aps_era_ECMWF_download(era_workdir,lonmin,lonmax,latmin,latmax,UTC,SLC1,SLC2)
|