Hi,
I have a netcdf file (time, lon, lat) on which I’d like to run these operations :
-
For each latitude of a subdomain, run a FFT2 : f(time,longitude) => ˆf(frequency,wavenumber) with
- the total spectrum
- the symmetrical part of the spectrum
- the anti-symmetrical part of the spectrum
-
Filter the summed total spectrum from (2) a little
-
Divide (1) by (2) and plot
For (1)
I think grdfft
does the job without specific options, so I assume I just need to run a loop such as
for lat in {lat1..lat2}; do
gmt grdcut file.nc -Rlon1/lon2/$lat/$lat -JM10c -Gsubset_$lat.nc
gmt grdfft subset_$lat.nc -Goutput_$lat.nc
done
Then for
- the total spectrum :
=> mean( Z{lat1:lat2} )
gmt grdmath output_*.nc -S MEAN = total.nc
- the anti-symmetric part :
=> mean( [Z{lat1:0} - Z{0:lat2}]/2 )
for lat in {lat1...0}; do
gmt grdmath output_$lat.nc ouput_$($lat+xx).nc SUB 2 DIV = anti_$lat.nc
done
gmt grdmath anti_*.nc -S MEAN = anti.nc
- the symmetric part :
=> mean( [Z{lat1:0} + Z{0:lat2}]/2 )
for lat in {lat1...0}; do
gmt grdmath output_$lat.nc ouput_$($lat+xx).nc ADD 2 DIV = symm_$lat.nc
done
gmt grdmath symm_*.nc -S MEAN = symm.nc
Is this even correctly translated to grdmath ?
For (2)
a simple use of grdfilter with -Fb
I assume, but which option corresponds best to a 1-2-1 filter ?
For (3)
And finally
gmt grdmath anti.nc filtered.nc DIV = anti_norm.nc
gmt grdmath symm.nc filtered.nc DIV = symm_norm.nc
gmt begin diagram png
gmt subplot 1x3 -R-15/15/0.001/0.5 -JX10c/10c+l
gmt subplot set
gmt grdimage filtered.nc
gmt subplot set
gmt grdimage anti_norm.nc
gmt subplot set
gmt grdimage symm_norm.nc
gmt subplot end
gmt end show
Do you have any input on steps (1) and (2) ?
Thanks,
[EDIT] : I forgot to mention an important part … the process is supposed to use a moving window of XX days for the FFT…