Commit 29548ad8 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merged from master

parents aac39c07 13e24f0d
#!/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,11 +86,20 @@ mm1=int(month)
dd1=int(day)
tt1=int(t1)
tt2=int(t2)
dT=3600*(tt1+(tt2-tt1)/2)
#dT=3600*(tt1+(tt2-tt1)/2)
#tte = dt.datetime(yy1, mm1, dd1, tt2)
ttt = dt.datetime(yy1, mm1, dd1, tt1)
tte = ttt + DT
#dT=3600*(tt2-tt1)
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)
......@@ -74,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])
......@@ -86,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]
......@@ -152,6 +187,7 @@ def make_netcdf(filename):
time.standard_name = 'time'
time.long_name = 'reference time for the measurements'
time.units = 'seconds since 2010-01-01 00:00:00'
time.bounds = 'time_bnds'
time[:] = timeA
time_bnds = ncf.createVariable('time_bnds','i4',('time_bnds',))
......@@ -161,31 +197,40 @@ def make_netcdf(filename):
time_bnds[:] = time_bndsA
#dBT = ncf.createVariable('dBT','f4',('lat','lon'),fill_value=9.96921e+36)
dBT = ncf.createVariable('dBT','f4',('lat','lon'))
dBT = ncf.createVariable('dBT','f4',('time', 'lat','lon'))
dBT.standard_name = 'dBT'
dBT.long_name = 'Brightness temperature difference (1231.5-1097.25 cm-1)'
dBT.units = 'K'
dBT[:,:] = comp1
dBT.ancillary_variables = 'flag'
dBT[0,:,:] = comp1
dBT_SeviriBand = ncf.createVariable('dBT_SeviriBand','f4',('lat','lon'))
dBT_SeviriBand = ncf.createVariable('dBT_SeviriBand','f4',('time','lat','lon'))
dBT_SeviriBand.standard_name = 'dBT'
dBT_SeviriBand.long_name = 'Brightness temperature difference (10-12 cm-1)'
dBT_SeviriBand.long_name = 'Brightness temperature difference (10-12 um)'
dBT_SeviriBand.units = 'K'
dBT_SeviriBand[:,:] = comp2
dBT_SeviriBand.ancillary_variables = 'flag'
dBT_SeviriBand[0,:,:] = comp2
dBT_flag = ncf.createVariable('dBT_flag','i4',('lat','lon'),fill_value=-9999)
dBT_flag = ncf.createVariable('dBT_flag','i4',('time','lat','lon'),fill_value=-9999)
dBT_flag.standard_name = 'status_flag'
dBT_flag.long_name = 'Ash detection flag'
dBT_flag.flag_values = 0, 1, 2
dBT_flag.flag_meanings = 'No_ash@dBT.lt.+1.5 Uncertain Ash@dBT.gt.+2.5'
dBT_flag[:,:] = comp3
dBT_flag[0,:,:] = comp3
dBT_SeviriBand_flag = ncf.createVariable('dBT_SeviriBand_flag','i4',('lat','lon'),fill_value=-9999)
dBT_SeviriBand_flag = ncf.createVariable('dBT_SeviriBand_flag','i4',('time','lat','lon'),fill_value=-9999)
dBT_SeviriBand_flag.standard_name = 'status_flag'
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[:,:] = 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()
......@@ -231,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]
......@@ -248,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)
......@@ -262,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
......@@ -278,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='/home/toves/Askevarsel/Latlon/'
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
\ No newline at end of file
#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, 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
sensX=sensX/2.0
sensY=sensY/2.0
! 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
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 (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
avg3=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
! TODO : these are not used?
avg1=sum(res_comp1(1:II))/II
avg2=sum(res_comp2(1:II))/II
avg3=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