Commit ea3ffb1d authored by Espen Sollum's avatar Espen Sollum
Browse files

merged all developments from satdata

parent 86e29f8a
#!/usr/bin/env python3.7
import matplotlib
matplotlib.use('Agg')
import numpy as np
......@@ -7,6 +8,7 @@ import scipy.constants as constants
#from Scientific.IO.NetCDF import NetCDFFile as Dataset
from netCDF4 import Dataset, num2date, date2num
import datetime as dt
from datetime import timedelta
import os
import sys
import matplotlib.pyplot as plt
......@@ -15,12 +17,32 @@ from matplotlib.ticker import MultipleLocator, FormatStrFormatter
import math
import pdb
# The aggregation time interval is set in makeCombined and should be read from the
# files in /Combined, along with the time of the last file (to be used for upper
# time_bnds)
DT = timedelta(hours=6)
latmin=45
latmax=85
lonmin=-50
lonmax=50
gridresX=0.25
gridresY=0.25
# gridresX=0.25
# gridresY=0.25
# gridresX=0.1
# gridresY=0.1
# Read grid resolution from './fortran_info_res.txt'
try:
with open('./fortran_info_res.txt', 'r') as f:
txt=f.readlines()
except:
print("ERROR: file not found:\n", os.getcwd() + '/fortran_info_res.txt')
sys.exit(0)
gridresX=float(txt[0])
gridresY=float(txt[1])
#Dimensions
dimLon = int((lonmax-lonmin)/gridresX+1)
......@@ -37,16 +59,18 @@ LAT =[]
LON =[]
cc1 =[]
cc2 =[]
cc3 =[]
with open(filename) as fid:
header=fid.readline()
print(header)
# print(header)
for line in fid:
data=line.split()
LAT.append(data[0])
LON.append(data[1])
cc1.append(data[2])
cc2.append(data[3])
cc2.append(data[3])
cc3.append(data[4])
# Find date from header
year=header[0:4]
......@@ -62,17 +86,20 @@ mm1=int(month)
dd1=int(day)
tt1=int(t1)
tt2=int(t2)
#dT=3600*(tt1+(tt2-tt1)/2)
tte=dt.datetime(yy1, mm1, dd1, tt2)
#tte = dt.datetime(yy1, mm1, dd1, tt2)
ttt = dt.datetime(yy1, mm1, dd1, tt1)
tte = ttt + DT
#dT=3600*(tt2-tt1)
dT = (tte-ttt).seconds
dT = (tte-ttt).total_seconds()
time0 = float((ttt - dt.datetime(2010, 1, 1)).total_seconds())
timeA = time0 + dT/2.0
timeA = time0 # + dT/2.0
time_bndsA = [time0, time0+dT]
#print(timeA)
#print(time_bndsA)
......@@ -80,6 +107,7 @@ time_bndsA = [time0, time0+dT]
#NB! Alternative, use pandas
comp1 = np.zeros([dimLat,dimLon])
comp2 = np.zeros([dimLat,dimLon])
compML = np.zeros([dimLat,dimLon])
comp3 = np.zeros([dimLat,dimLon])
comp4 = np.zeros([dimLat,dimLon])
llat = np.zeros([dimLat,dimLon])
......@@ -92,6 +120,7 @@ for i in range(0,dimLon,1):
ii=ii+1
comp1[j,i]=cc1[ii]
comp2[j,i]=cc2[ii]
compML[j,i]=cc3[ii]
llat[j,i]=LAT[ii]
llon[j,i]=LON[ii]
......@@ -194,7 +223,14 @@ def make_netcdf(filename):
dBT_SeviriBand_flag.long_name = 'Ash detection flag (Seviri band)'
dBT_SeviriBand_flag.flag_values = 0, 1, 2
dBT_SeviriBand_flag.flag_meanings = 'No_ash@dBT.gt.-0.2 Uncertain Ash@dBT.lt.-0.5)'
dBT_SeviriBand_flag[0,:,:] = comp4
dBT_SeviriBand_flag[0,:,:] = comp4
AshML = ncf.createVariable('MassLoadingSEVIRI','f4',('time', 'lat','lon'))
AshML.standard_name = 'MassLoadingSEVIRI'
AshML.long_name = 'Ash Mass loading, SEVIRI method'
AshML.units = 'g m^-2'
AshML[0,:,:] = compML
ncf.close()
......@@ -240,8 +276,14 @@ def PlotLatLon(type,plotfile):
ylab='Flag (0=no_ash,1=uncertain,2=ash)'
plot_title=('IASI flag (Seviri band) '+ header)
comp=comp4
if (type<3):
elif(type==5):
# Plot mass loading
lim=[0, 10]
ylab='Ash column mass'
plot_title=('Ash column mass (SEVIRI) '+ header)
comp=compML
if (type<3 or type==5):
mycmap=plt.cm.jet
else:
levels = [0, 1, 2, 3]
......@@ -257,7 +299,7 @@ def PlotLatLon(type,plotfile):
plt.title(plot_title)
h=plt.colorbar(label=ylab)
if (type>2):
if (type>2 and type != 5):
#h.ax.locator_params(nbins=3)
h.set_ticks([0.33, 1, 1.66])
h.ax.set_yticklabels([0,1,2], fontsize=12)
......@@ -271,7 +313,7 @@ def PlotLatLon(type,plotfile):
##################################################
#Create netcdf file
filename_nc='IASI_grid_'+year+month+day+'_'+t1+'-'+t2+'.nc'
print(filename_nc)
#print(filename_nc)
make_netcdf(filename_nc)
#Create plot and filename
......@@ -287,3 +329,6 @@ PlotLatLon(3,outfilep)
outfilep='IASI_Sflag_'+year+month+day+'_'+t1+'-'+t2+'.png'
PlotLatLon(4,outfilep)
outfilep='IASI_SMass_'+year+month+day+'_'+t1+'-'+t2+'.png'
PlotLatLon(5,outfilep)
#!/usr/bin/env python3.7
"""
Extract data from a netcdf variable
"""
import sys
import numpy as np
from netCDF4 import Dataset
file_name = sys.argv[1]
var_name = sys.argv[2]
data = Dataset(file_name, 'r')
var = data.variables[var_name][:]
np.savetxt(var_name + '.txt', var,fmt="%f")
......@@ -3,8 +3,20 @@
### Run script (bash, fortran and python):
### ./regrid.sh Combined_IASI_20100417_1200-2400.nc
path_in=./
file=$1
# path_in=./
# file=$1
# Exit if not 2 args given (path file)
if [[ $# -ne 2 ]] ; then
echo 'usage: regrid.sh path_in file'
exit 0
fi
path_in=$1
file=$2
# Go to work directory
cd ${path_in}
year=${file:14:4}
month=${file:18:2}
......@@ -26,15 +38,27 @@ time=${d1}:00-${d2}:00
#ncks -s "%e\n" -C -H -v lon Combined_IASI_20191128_1800-2400.nc > lon.txt
#ncks -s "%e\n" -C -H -v dBT Combined_IASI_20191128_1800-2400.nc > dBT.txt
echo $path_in$file
#echo $path_in$file
# ncks -s "%e\n" -C -H -v lat $path_in$file > lat.txt
# ncks -s "%e\n" -C -H -v lon $path_in$file > lon.txt
# ncks -s "%e\n" -C -H -v dBT $path_in$file > dBT.txt
# ncks -s "%e\n" -C -H -v SEVIRI_dBT $path_in$file > SEVIRI_dBT.txt
# Can not use ncks on betty
python3.7 extract_nc.py $path_in$file lat # > lat.txt
python3.7 extract_nc.py $path_in$file lon # > lon.txt
python3.7 extract_nc.py $path_in$file dBT # > dBT.txt
python3.7 extract_nc.py $path_in$file SEVIRI_dBT # > SEVIRI_dBT.txt
python3.7 extract_nc.py $path_in$file MassLoading # > MassLoadingSEVIRI.txt
mv MassLoading.txt MassLoadingSEVIRI.txt
# TODO: rename current 'MassLoading' for IASI/SEVIRI mass
ncks -s "%e\n" -C -H -v lat $path_in$file > lat.txt
ncks -s "%e\n" -C -H -v lon $path_in$file > lon.txt
ncks -s "%e\n" -C -H -v dBT $path_in$file > dBT.txt
ncks -s "%e\n" -C -H -v SEVIRI_dBT $path_in$file > SEVIRI_dBT.txt
GridresolutionX=0.25
GridresolutionY=0.25
# GridresolutionX=0.1
# GridresolutionY=0.1
sensX=2.0
sensY=0.5
latmin=45
......@@ -42,7 +66,7 @@ latmax=85
lonmin=-50
lonmax=50
rm fortran_info_res.txt
rm -f fortran_info_res.txt
echo $GridresolutionX > fortran_info_res.txt
echo $GridresolutionY >> fortran_info_res.txt
echo $sensX >> fortran_info_res.txt
......@@ -59,9 +83,10 @@ lonmax=50
echo lon.txt >> fortran_info_res.txt
echo dBT.txt >> fortran_info_res.txt
echo SEVIRI_dBT.txt >> fortran_info_res.txt
echo MassLoadingSEVIRI.txt >> fortran_info_res.txt
./regrid_IASI
python3 MakeIASI_netcdf_and_Plot.py
python3.7 MakeIASI_netcdf_and_Plot.py
rm lat.txt lon.txt dBT.txt SEVIRI_dBT.txt
#rm -f lat.txt lon.txt dBT.txt SEVIRI_dBT.txt
Program regrid_IASI
Implicit none
Integer I,J,ILINE,IO,II, INDEX1, INDEX2
real :: mm,dd,avg1, avg2
real :: dx, dy,sensX, sensXorg, sensY, min_err, GridresX, GridresY, rlon,rlat
character*50 :: path
character FILE1*20,FILE2*20,FILE3*20,FILE4*20
character FILEgrid*17,yyt*4, mmt*2, ddt*2, time*11
real:: lat(100000),lon(100000),comp1(100000), comp2(100000)
real:: res_lat(8000),res_lon(8000),res_err(8000),res_comp1(8000)
real:: res_comp2(8000)
integer :: ilat, ilon, qa, qa_lim
integer :: nx, ny
real :: latmin,latmax,lonmin,lonmax
character(len=10) :: fmt ! format descriptor
program regrid_IASI
Implicit none
Integer I,J,ILINE,IO,II, INDEX1, INDEX2
real :: mm,dd,avg1, avg2, avg3
real :: dx, dy,sensX, sensXorg, sensY, min_err, GridresX, GridresY, rlon,rlat
character*50 :: path
character FILE1*20,FILE2*20,FILE3*20,FILE4*20,FILE5*30
character FILEgrid*17,yyt*4, mmt*2, ddt*2, time*11
real:: lat(100000),lon(100000),comp1(100000), comp2(100000), comp3(100000)
real:: res_lat(8000),res_lon(8000),res_err(8000),res_comp1(8000)
real:: res_comp2(8000), res_comp3(8000)
integer :: ilat, ilon, qa, qa_lim
integer :: nx, ny
real :: latmin,latmax,lonmin,lonmax
logical :: ldum=.false.
character(len=10) :: fmt ! format descriptor
!ncks -s "%e\n" -C -H -v lat Combined_IASI_20191128_1800-2400.nc > lat.txt
!ncks -s "%e\n" -C -H -v lon Combined_IASI_20191128_1800-2400.nc > lon.txt
!ncks -s "%e\n" -C -H -v dBT Combined_IASI_20191128_1800-2400.nc > dBT.txt
FILEgrid='YYYYMMDD_grid.out'
open(11,File='fortran_info_res.txt', status='old')
IO=0
read(11,*,END=50) GridresX
read(11,*,END=50) GridresY
read(11,*,END=50) sensX
read(11,*,END=50) sensY
read(11,*,END=50) latmin
read(11,*,END=50) latmax
read(11,*,END=50) lonmin
read(11,*,END=50) lonmax
read(11,'(a4)',END=50) yyt
read(11,'(a2)',END=50) mmt
read(11,'(a2)',END=50) ddt
read(11,'(a11)',END=50) time
read(11,*,END=50) IO,I,I,FILE1
read(11,*,END=50) FILE2
read(11,*,END=50) FILE3
read(11,*,END=50) FILE4
FILEgrid='YYYYMMDD_grid.out'
open(11,File='fortran_info_res.txt', status='old')
IO=0
read(11,*,END=50) GridresX
read(11,*,END=50) GridresY
read(11,*,END=50) sensX
read(11,*,END=50) sensY
read(11,*,END=50) latmin
read(11,*,END=50) latmax
read(11,*,END=50) lonmin
read(11,*,END=50) lonmax
read(11,'(a4)',END=50) yyt
read(11,'(a2)',END=50) mmt
read(11,'(a2)',END=50) ddt
read(11,'(a11)',END=50) time
read(11,*,END=50) IO,I,I,FILE1
read(11,*,END=50) FILE2
read(11,*,END=50) FILE3
read(11,*,END=50) FILE4
read(11,*,END=50) FILE5
!iline=IO/3 - 2
iline=IO -2
iline=IO -2
! Input files have only one data point -> treat as dummy file
if(IO == 1) ldum=.true.
! eso: scale sensX with latitude (keep roughly original values at latmin)
sensXorg=sensX/2.0*abs(cos(latmin*3.14/180.))
sensY=sensY/2.0
sensXorg=sensX/2.0*abs(cos(latmin*3.14/180.))
sensY=sensY/2.0
50 continue
close(11)
close(11)
IO=1
IO=1
!satfile=trim(path)//trim(FILE1)
write(*,*) FILE1
write(*,*) FILE2
write(*,*) FILE3
write(*,*) FILE4
write(*,*) 'iline: ',iline
lat(:)=0.0
lon(:)=0.0
comp1(:)=0.0
! write(*,*) FILE1
! write(*,*) FILE2
! write(*,*) FILE3
! write(*,*) FILE4
! write(*,*) 'iline: ',iline
lat(:)=0.0
lon(:)=0.0
comp1(:)=0.0
comp2(:)=0.0
comp3(:)=0.0
open(11,File=trim(FILE1), status='old')
do j=1,iline
read(11,*) lat(j)
read(11,*) lat(j)
end do
close(11)
open(11,File=trim(FILE2), status='old')
do j=1,iline
read(11,*) lon(j)
if (lon(j).ge.180.0) lon(j)=lon(j)-360.0
end do
read(11,*) lon(j)
if (lon(j).ge.180.0) lon(j)=lon(j)-360.0
end do
close(11)
open(11,File=trim(FILE3), status='old')
do j=1,iline
read(11,*) comp1(j)
end do
read(11,*) comp1(j)
end do
close(11)
open(11,File=trim(FILE4), status='old')
do j=1,iline
read(11,*) comp2(j)
end do
read(11,*) comp2(j)
end do
close(11)
open(11,File=trim(FILE5), status='old')
do j=1,iline
read(11,*) comp3(j)
end do
close(11)
open(12,File=FILEgrid, status='unknown')
write(12,'(a4,a1,a2,a1,a2,1x,a11)')yyt,'-',mmt,'-',ddt,time
open(12,File=FILEgrid, status='unknown')
write(12,'(a4,a1,a2,a1,a2,1x,a11)')yyt,'-',mmt,'-',ddt,time
!############################
nx = nint((lonmax - lonmin)/GridresX)
ny = nint((latmax - latmin)/GridresY)
do ilon=0,nx
do ilat=0,ny
rlon = lonmin + ilon*gridresX
rlat = latmin + ilat*gridresY
II=0
res_lat(:)=0.0
res_lon(:)=0.0
res_err(:)=999.0
res_comp1(:)=0.0
INDEX1=0
INDEX2=0
avg1=0.0
avg2=0.0
do j=1,iline
if (abs(lat(J)) < 89.) sensX = sensXorg/abs(cos(lat(J)*3.14/180.))
if (lat(J).ge.(rlat-sensY) .and. lat(J).le.(rlat+sensY) .and. &
lon(J).ge.(rlon-sensX) .and. lon(J).le.(rlon+sensx)) then
dx=(lat(J)-rlat)*(lat(J)-rlat)
dy=(lon(J)-rlon)*(lon(J)-rlon)
II=II+1
res_lat(II)=lat(J)
res_lon(II)=lon(J)
res_err(II)=sqrt(dx+dy)
res_comp1(II)=comp1(J)
res_comp2(II)=comp2(J)
end if
end do
if (II.gt.0) then
avg1=SUM(res_comp1(1:II))/II
avg2=SUM(res_comp2(1:II))/II
min_err=minval(res_err) !Closest grid
INDEX1=MINLOC(res_err,DIM=1) !Index to corresponding minumum
!INDEX2=MINLOC(res_err,DIM=1,MASK=res_comp1.GT.0) !Index to grid with comp1 > 0
write(12,'(4f12.2,I6)') rlat,rlon, res_comp1(INDEX1),res_comp2(INDEX1), II
else
!write(12,'(3f12.2,I6)') rlat,rlon, -99.9, 0
write(12,'(2f12.2,2E13.6,I6)') rlat,rlon, 9.96921e36, 9.96921e36,0
end if
end do !lat
end do !lon
end
nx = nint((lonmax - lonmin)/GridresX)
ny = nint((latmax - latmin)/GridresY)
if (ldum.eqv..false.) then
do ilon=0,nx
do ilat=0,ny
rlon = lonmin + ilon*gridresX
rlat = latmin + ilat*gridresY
II=0
res_lat(:)=0.0
res_lon(:)=0.0
res_err(:)=999.0
res_comp1(:)=0.0
INDEX1=0
INDEX2=0
avg1=0.0
avg2=0.0
do j=1,iline
if (abs(lat(J)) < 89.) sensX = sensXorg/abs(cos(lat(J)*3.14/180.))
if (lat(J).ge.(rlat-sensY) .and. lat(J).le.(rlat+sensY) .and. &
lon(J).ge.(rlon-sensX) .and. lon(J).le.(rlon+sensx)) then
dx=(lat(J)-rlat)*(lat(J)-rlat)
dy=(lon(J)-rlon)*(lon(J)-rlon)
II=II+1
res_lat(II)=lat(J)
res_lon(II)=lon(J)
res_err(II)=sqrt(dx+dy)
res_comp1(II)=comp1(J)
res_comp2(II)=comp2(J)
res_comp3(II)=comp3(J)
end if
end do
if (II.gt.0) then
avg1=sum(res_comp1(1:II))/II
avg2=sum(res_comp2(1:II))/II
avg2=sum(res_comp3(1:II))/II
min_err=minval(res_err) !Closest grid
INDEX1=MINLOC(res_err,DIM=1) !Index to corresponding minumum
!INDEX2=MINLOC(res_err,DIM=1,MASK=res_comp1.GT.0) !Index to grid with comp1 > 0
write(12,'(5f12.2,I6)') rlat,rlon, res_comp1(INDEX1),res_comp2(INDEX1), res_comp3(INDEX1), II
else
!write(12,'(3f12.2,I6)') rlat,rlon, -99.9, 0
write(12,'(2f12.2,3E13.6,I6)') rlat,rlon, 9.96921e36, 9.96921e36, 9.96921e36, 0
end if
end do !lat
end do !lon
else
! Dummy file, write fill values everywhere
do ilon=0,nx
do ilat=0,ny
write(12,'(2f12.2,3E13.6,I6)') rlat,rlon, 9.96921e36, 9.96921e36, 9.96921e36, 0
end do
end do
end if
close(12)
end program regrid_IASI
#!/usr/bin/env python3.7
"""
Run IASI regridding scripts
"""
import sys
import os
import subprocess
import argparse
from datetime import timedelta, datetime
import time
import re
import shutil
import atexit
def at_exit(file_name):
""" Delete lock file on unexpected exit """
try:
os.unlink(file_name)
except Exception as any_err:
pass
def get_src_files(src_dir, tmp_dir, days_tp):
"""
Scan the source directory for all available IASI composite files to be regridded.
Copy eligble files to tmp directory.
Filter by days_tp (days to process).
"""
list_files_all = [f for f in os.listdir(src_dir) if re.match(r'Combined_IASI_.*\.nc', f)]
# Files to process
file_list = []
t_now = datetime.utcnow()
for fn in list_files_all:
t_file = str_to_dt(fn)
if (t_now - t_file).days <= days_tp:
file_list.append(fn)
# Copy to tmp directory
for f in file_list:
shutil.copyfile(os.path.join(src_dir, f), os.path.join(tmp_dir, f))
return file_list
def str_to_dt(fn):