Created
March 23, 2023 13:12
-
-
Save Gro-Tsen/c40b63c61d895aed0ce5f9168027debb to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #### Plot of cumulative precipitations in France | |
| ## Input data: /tmp/precs.csv is expected to contain a list of daily | |
| ## precipitations in France (first column = date (unused), second | |
| ## column = precipiptations in mm/day). It can be produced from the | |
| ## data source described in <URL: | |
| ## https://twitter.com/gro_tsen/status/1620435663844970497 > by the | |
| ## following command line: | |
| # perl -ne 'next if /^\#/; die unless /^([0-9]{8})\s+(\-?[0-9]+\.[0-9]+(?:[Ee](?:[+-]?[0-9]+))?)\s*$/; printf "%s,%.4f\n", $1,($2+0)' /data/meteo/climexp/iera5_prcp_daily_eu_France_metropolitan.dat > /tmp/precs.csv | |
| ### Read data | |
| import csv | |
| data = list(csv.reader(open("/tmp/precs.csv","r"))) | |
| data = list(map(lambda x: [x[0],float(x[1])],data)) | |
| datamean = mean([v[1] for v in data]) | |
| datayears = floor(len(data)/365.25) | |
| # Add 10 years of constant data at the start that will be removed after convolution: | |
| fakedatacnt = ceil((datayears+10)/4)*1461 - len(data) | |
| extdata = [datamean]*fakedatacnt + [v[1] for v in data] | |
| extdatayears = len(extdata)/1461*4 | |
| ### Convolve data with exponentials (by multiplying their Fourier transforms) | |
| fft = FastFourierTransform(len(extdata)) | |
| fftconvol = FastFourierTransform(len(extdata)) | |
| fftconvolq = FastFourierTransform(len(extdata)) | |
| for i in range(len(extdata)): | |
| fft[i] = (extdata[i], 0) | |
| fftconvol[i] = (exp(-i/365.25), 0) | |
| fftconvolq[i] = (exp(-i/(365.25/4)), 0) | |
| fft.forward_transform() | |
| fftconvol.forward_transform() | |
| fftconvolq.forward_transform() | |
| fftout = FastFourierTransform(len(extdata)) | |
| fftoutq = FastFourierTransform(len(extdata)) | |
| for i in range(len(extdata)): | |
| before = CC(fft[i][0], fft[i][1]) | |
| coeft = (CC(fftconvol[i][0], fftconvol[i][1]))*365.25/fftconvol[0][0] | |
| # Should be roughly equal to: | |
| # 365.25 * (1/(1+2*I*pi*i/extdatayears) + 1/(1+2*I*pi*(i-len(extdata))/extdatayears)) | |
| after = before*coeft | |
| fftout[i]=(real(after), imag(after)) | |
| coeftq = (CC(fftconvolq[i][0], fftconvolq[i][1]))*(365.25/4)/fftconvolq[0][0] | |
| afterq = before*coeftq | |
| fftoutq[i]=(real(afterq), imag(afterq)) | |
| fftout.inverse_transform() | |
| fftoutq.inverse_transform() | |
| # Convolved data: | |
| convold = [(data[i][0], fftout[i+fakedatacnt][0]) for i in range(len(data))] | |
| convoldq = [(data[i][0], fftoutq[i+fakedatacnt][0]) for i in range(len(data))] | |
| ### Perform long-term regression on data | |
| import scipy.stats | |
| regr = scipy.stats.linregress([(i/365.25,data[i][1]) for i in range(len(data))]) | |
| ### Plot the result | |
| g = plot((regr.slope*(x-1950-0.5)+regr.intercept)*365.25, (x,1951,1950+len(data)/365.25), color="red") + line([(i/365.25+1950, convold[i][1]) for i in range(365,len(data))], title="Précipitation cumulée en France (sur 1 an à amort. exp., en mm/an)") | |
| g.show(dpi=192) | |
| gq = plot((regr.slope*(x-1950-0.125)+regr.intercept)*365.25/4, (x,1951,1950+len(data)/365.25), color="red") + line([(i/365.25+1950, convoldq[i][1]) for i in range(365,len(data))], title="Précipitation cumulée en France (sur ¼ an à amort. exp., en mm/(¼an))", thickness=0.4) | |
| gq.show(dpi=192) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment