From 6e685760b9b45267ae9025558d32a8251befb159 Mon Sep 17 00:00:00 2001 From: Paco Lopez Dekker <f.lopezdekker@tudelft.nl> Date: Thu, 18 Apr 2024 15:00:32 +0200 Subject: [PATCH] Import of existing CoSAR code --- .DS_Store | Bin 0 -> 6148 bytes .gitignore | 2 + cosar/__init__.py | 16 ++ cosar/geosim1d.py | 316 +++++++++++++++++++++++++++++ cosar/geosim1d_test.py | 50 +++++ cosar/results.py | 54 +++++ cosar/simplesim.py | 445 +++++++++++++++++++++++++++++++++++++++++ cosar/viciprep.py | 67 +++++++ cosarsim.py | 40 ++++ setup.cfg | 27 +++ setup.py | 15 ++ 11 files changed, 1032 insertions(+) create mode 100644 .DS_Store create mode 100644 .gitignore create mode 100644 cosar/__init__.py create mode 100644 cosar/geosim1d.py create mode 100644 cosar/geosim1d_test.py create mode 100644 cosar/results.py create mode 100644 cosar/simplesim.py create mode 100644 cosar/viciprep.py create mode 100644 cosarsim.py create mode 100644 setup.cfg create mode 100644 setup.py diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..abe499281aa11f982c6ff7a17d65ffecdf6ccaae GIT binary patch literal 6148 zcmeH~I|>3p42BaQAlO)1PU8W*!654iynwG#SWxWe=>B;ixLS+IA4vY0Orq>p>}*6t zH}~^eWFaCmxT&lx3{0_~%1L^;$w98?+v#%b`-@kNves(g_>I@|Jf@HU36KB@kN^q% z5COZlVY7KCBMFcI2|Nkd`=P*1YibMiR|kTR0MH52Zdm&)0WFq**3=e?3{0aH8m;PM zh~>Q<TC%RDw$Nx7&EZ4ypVg)qm`1y3!3L()g@FV}U_f9T`-R>ATllB>f6&4$36Q{_ z5ztAu>o$0(yj$NM&+_}s+PcA^UXJkc5rB<d#Y?yw&WkOeHMNBz1LKE)V_+bGuM&6w DS@RQ) literal 0 HcmV?d00001 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fd20fdd --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ + +*.pyc diff --git a/cosar/__init__.py b/cosar/__init__.py new file mode 100644 index 0000000..f5f9e97 --- /dev/null +++ b/cosar/__init__.py @@ -0,0 +1,16 @@ +""" +=================================== +Input and output (:mod:`pycosar`) +=================================== + +.. currentmodule:: pycosar + + +""" + +from simplesim import * +from results import * +from geosim1d import * +#from rat import * + +__all__ = ['simplesim','results'] \ No newline at end of file diff --git a/cosar/geosim1d.py b/cosar/geosim1d.py new file mode 100644 index 0000000..d205180 --- /dev/null +++ b/cosar/geosim1d.py @@ -0,0 +1,316 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Dec 18 16:13:58 2013 + +@author: lope_fr +""" + + +import numpy as np +import matplotlib +import matplotlib.pyplot as plt +from scipy import constants +from drama import utils +# csl stands fro cosarlib +import cosar.simplesim as csl +import matplotlib +matplotlib.rcParams['ps.fonttype'] = 42 + +def geosim1d(cosar, surface, dx, n_rg=1, PRF=0, Tint=300.0, hamming=True, + n_lags=1, surfaceseed=None, compute_IRF=False, pos_pt=[]): + """ + Simple 1D simulation using the geometry defined in cosar + @type cosar csl.CoSARGeometry + """ + #Init cosar if nothing given. The secquence also forces the + #assert isinstance(cosar, csl.CoSARGeometry) + PRF_min = 2*cosar.dDoppler_bandwidth() + if PRF == 0: + PRF = PRF_min + else: + if PRF < PRF_min: + print("PRF given to low") + PRF = PRF_min + #Check that the azimuth resolution is finer than the minimum lenth scale + if dx > surface.L_min()/2: + dx = surface.L_min()/2 + n_az = int(np.round(cosar.L_illuminated()/dx)) + n_t = int(Tint * PRF) + msg = "Number of azimuth grid samples: " + repr(n_az) + print(msg) + msg = "Number of pulses: " + repr(n_t) + print(msg) + k0 = 2*np.pi*cosar.f0/constants.c + #Seed for repeatibility of the surface + np.random.seed(surfaceseed) + nrcs, dop, h, az = surface.scene(cosar.L_illuminated(), dx, cosar.f0) + if compute_IRF: + nrcs[:] = 0 + if np.size(pos_pt) == 0: + pos_pt = [0] + #nrcs[n_az/2 + np.round(1e3/dx)] = 1 + dop[:] = 0 + h[:] = 0 + for pos_this_pt in pos_pt: + pos_this_pt_s = np.round(pos_this_pt/dx) + n_az/2 + nrcs[pos_this_pt_s] = 1 + #Reseed to randomize realization + np.random.seed() + #Time vector + t_v = np.linspace(-Tint / 2, Tint / 2, int(Tint * PRF)) + #Doppler phase shift for each point on the ground + dop_phasor = np.exp(2j * np.pi * dop / PRF) + #The expected amplitude takes into account the grid size and the fact + #that we are working with circular Gaussian signals + expected_amp = np.sqrt(nrcs * dx / 2) + #Now we calculate the trajectory of the two spacecraft assuming linear + #motions + pos1 = np.zeros([n_t, 3]) + pos2 = np.zeros([n_t, 3]) + pos1[:, 1] = -np.sin(cosar.theta_i) * cosar.R0() + pos2[:, 1] = -np.sin(cosar.theta_i - cosar.dtheta_i()) * cosar.R0() + pos1[:, 2] = np.cos(cosar.theta_i) * cosar.R0() + pos2[:, 2] = np.cos(cosar.theta_i - cosar.dtheta_i()) * cosar.R0() + for dim_ind in range(3): + pos1[:, dim_ind] = pos1[:, dim_ind] + cosar.dv[dim_ind] / 2.0 * t_v + pos2[:, dim_ind] = pos2[:, dim_ind] - cosar.dv[dim_ind] / 2.0 * t_v + #Correlation coefficcient + # alpha**(tau*PRF) = 1/e + # tau*PRF * log(alpha)=-1 alpha = exp(-1/(tau*PRF) + alpha = 1.0 + if surface.tau > 0: + alpha = np.exp(-1.0/(surface.tau*PRF)) + msg = "Pulse to pulse coherence: " + repr(alpha) + print(msg) + alpha_end = np.exp(-1.0*t_v.size/(surface.tau*PRF)) + msg = "Start to end coherence: " + repr(alpha_end) + print(msg) + if n_rg > 1: + uscat = np.random.randn(n_az, n_rg) + 1j*np.random.randn(n_az, n_rg) + else: + uscat = np.sqrt(2) * np.exp(1j*2*np.pi*np.random.rand(n_az, n_rg)) + #Now main loop + for t_ind in range(0, t_v.size): + #Generate random complex scatterers + uscat = (alpha * uscat * dop_phasor.reshape(n_az, 1) + + np.sqrt(1-alpha**2) * (np.random.randn(n_az, n_rg) + + 1j * np.random.randn(n_az, n_rg))) + #Multipliy them with a fixed amplitude + scat = uscat * expected_amp.reshape(n_az, 1) + #Calculate ranges from targets to both radars + r1 = ((np.sqrt((pos1[t_ind, 1])**2 + + (az - pos1[t_ind, 0])**2 + + (pos1[t_ind, 2] - h)**2)).reshape(n_az, 1)) + r2 = ((np.sqrt((pos2[t_ind, 1])**2 + + (az - pos2[t_ind, 0])**2 + + (pos2[t_ind, 2] - h)**2)).reshape(n_az, 1)) + #Include 2-way phase in field + s1 = scat * np.exp(-2j*k0*r1) + s2 = scat * np.exp(-2j*k0*r2) + # Compute range profiles by integrating in azimuth. + # This assumed that there is no range migration + rx1_ = s1.sum(axis=0) + rx2_ = s2.sum(axis=0) + # Save in a temporal stack + if t_ind == 0: + rx1 = rx1_.copy() + rx2 = rx2_.copy() + else: + rx1 = np.vstack((rx1, rx1_)) + rx2 = np.vstack((rx2, rx2_)) + #Now we have a range-slow time stack of CoSAR data + n_az_out = n_az + proc_img = np.zeros([2 * n_lags + 1, n_az_out], dtype=np.complex) + for lag_ind in range(0, 2*n_lags + 1): + lag = lag_ind - n_lags + proc_img[lag_ind, :], T_i_scaling = \ + csl.cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0, + han=hamming, dt=lag) + #Backprojection SAR focusing of one channel + sar_img, kk = csl.sar_proc_1d(rx1, pos1, az, k0, han=hamming) + #Amplitude scaling of processed image + C_p = 2 * cosar.dv[0] / (cosar.wavelength() * cosar.R0()) + amp_scaling = C_p / PRF + proc_img = amp_scaling * proc_img + sar_img = amp_scaling * sar_img + T_i_scaled = T_i_scaling * Tint + return (az, nrcs, dop, h, proc_img, n_lags, T_i_scaled, dx, n_rg, sar_img) + + +def geosim1d_montecarlo(cosar, surface, dx, n_rg=1, PRF=0, Tint=300.0, + hamming=True, n_lags=1, n_it=10, surfaceseed=1): + """ + A function to make several runs + """ + sim_res = geosim1d(cosar, surface, dx, n_rg=n_rg, PRF=PRF, Tint=Tint, + n_lags=n_lags, hamming=hamming, surfaceseed=surfaceseed) + #Unpack first simulation + az, nrcs, dop, h, proc_img, n_lags, T_i, dx, n_rg, sar_img = sim_res + dim_res = proc_img.shape + dim_mtc = (n_it,) + dim_res + proc_imgs = np.zeros(dim_mtc, dtype=np.complex) + #Insert output of first simulation in stack + proc_imgs[0, :, :] = proc_img + for it in range(1, n_it): + sim_res = geosim1d(cosar, surface, dx, n_rg=n_rg, PRF=PRF, Tint=Tint, + n_lags=n_lags, hamming=hamming, + surfaceseed=surfaceseed) + az, nrcs, dop, h, proc_img, n_lags, T_i, dx, n_rg, sar_img = sim_res + proc_imgs[it, :, :] = proc_img + + return (az, nrcs, dop, h, proc_imgs, n_lags, T_i, dx, n_rg, PRF) + + +def geosim1d_save(sim_results, fileroot): + """ + Function to save simulator results + """ + #Unpack simulation results + az, nrcs, dop, h, proc_imgs, n_lags, T_i, dx, n_rg, PRF = sim_results + filename = fileroot+'_res.npz' + np.savez(filename, az=az, nrcs=nrcs, dop=dop, h=h, proc_imgs=proc_imgs, + n_lags=n_lags, T_i=T_i, dx=dx, n_rg=n_rg, PRF=PRF) + + +def geosim1d_load(fileroot): + """ + Function to save simulator results + """ + filename = fileroot+'_res.npz' + res_d = np.load(filename) + res = (res_d['az'], res_d['nrcs'], res_d['dop'], res_d['h'], + res_d['proc_imgs'], res_d['n_lags'], res_d['T_i'], + res_d['dx'], res_d['n_rg'], res_d['PRF']) + return (res) + + +def geosim1d_analyze(cosar, surface, sim_results, dB=False, fileroot=False, + fontsize=16): + """ + Function to analyze output of geosim1d + """ + matplotlib.rcParams.update({'font.size': fontsize}) + #Unpack simulation results + az, nrcs, dop, h, proc_imgs, n_lags, T_i, dx, n_rg, PRF = sim_results + #Calculate Quality factor + az_res = 0.5 * cosar.wavelength() * cosar.R0() / (cosar.dv[0] * T_i) + #tau_ca + tau_ca_Doppler = (cosar.wavelength()*np.sqrt(2) / + (4 * np.pi * surface.vr_stdev)) + tau_ca = np.sqrt(1.0 / (1 / surface.tau**2 + 1 / tau_ca_Doppler**2)) + tau_ca_quant = np.max([tau_ca, 1 / PRF]) + Q_nrcs = csl.Q_w(nrcs, nrcs, dx, n_rg, T_i, tau_ca_quant, az_res) + #Expected standard deviation of the estimated nrcs + nrcs_err_t = nrcs / np.sqrt(Q_nrcs) + #Lag 1 expected value + dop_phase = 2 * np.pi * dop / PRF + dop_phasor = np.exp(2j * np.pi * dop / PRF) + #Height phasor + h_phase = cosar.kz() * h + h_phasor = np.exp(1j * h_phase) + corr_ati = 1.0 + if surface.tau > 0: + corr_ati = np.exp(-1.0 / (surface.tau * PRF)) + R_tau_lag1 = corr_ati * dop_phasor * h_phasor * nrcs + R_tau_lag1m = corr_ati * np.conj(dop_phasor) * h_phasor * nrcs + R_tau_lag0 = h_phasor * nrcs + if np.size(proc_imgs.shape) == 3: + proc_img = proc_imgs[0, :, :] + nrcs_err = proc_imgs[:, n_lags, :] - R_tau_lag0.reshape((1, nrcs.size)) + nrcs_err_mean = np.mean(nrcs_err, axis=0) + nrcs_err_std = np.std(nrcs_err, axis=0) + plt.figure() + ax = plt.subplot(1, 1, 1) + p_std, = plt.plot(az / 1e3, nrcs_err_std, 'b', + label=r'std($\delta R_{\tau}$)') + p_mean_r, = plt.plot(az / 1e3, np.real(nrcs_err_mean), 'c', + label=r'$\mathrm{Re}\{\overline{\delta R_{\tau}}\}$') + p_mean_i, = plt.plot(az / 1e3, np.imag(nrcs_err_mean), 'm', + label=r'$\mathrm{Im}\{\overline{\delta R_{\tau}}\}$') + plt.plot(az / 1e3, nrcs_err_t, 'b--') + plt.xlabel('Azimuth [km]', fontsize=fontsize) + plt.ylabel('NRCS errors', fontsize=fontsize) + ax.set_ylim([-0.05,0.1]) + handles, labels = ax.get_legend_handles_labels() + l1 = plt.legend([handles[0]], [labels[0]], loc=2) + l2 = plt.legend(handles[1:], labels[1:], loc=3, ncol=2) + plt.gca().add_artist(l1) + if type(fileroot) == str: + filename = fileroot+'_err_stats.eps' + plt.savefig(filename) + filename = fileroot+'_err_stats.png' + plt.savefig(filename) + plt.figure() + plt.plot(az/1e3, np.real(R_tau_lag1)) + for ind in range(1, proc_imgs.shape[0]): + plt.plot(az/1e3, np.real(proc_imgs[ind, n_lags+1, :])) + + else: + proc_img = proc_imgs + + #Estimate nrcs + nrcs_cosar = np.abs(proc_img[n_lags, :]) + #FIXME + #Show results + #az in km + az = az / 1e3 + #plt.figure() + #csl.plot_nrcs(az, nrcs, az, nrcs_cosar, dB=dB) + #Unwrap CoSAR using known true phase. Since this is not a phase unwrapping + #exercise, this is not cheating. + R_tau_lag0_CoSAR = proc_img[n_lags, :] + R_tau_lag0_CoSAR_a = np.abs(R_tau_lag0_CoSAR) + R_tau_lag0_CoSAR_p = (np.angle(R_tau_lag0_CoSAR * np.conj(R_tau_lag0)) + + h_phase) + R_tau_lag1_CoSAR = proc_img[n_lags+1, :] + R_tau_lag1_CoSAR_a = np.abs(R_tau_lag1_CoSAR) + R_tau_lag1_CoSAR_p = (np.angle(R_tau_lag1_CoSAR * np.conj(R_tau_lag1)) + + dop_phase + h_phase) + R_tau_lag1m_CoSAR = proc_img[n_lags-1, :] + R_tau_lag1m_CoSAR_a = np.abs(R_tau_lag1m_CoSAR) + R_tau_lag1m_CoSAR_p = (np.angle(R_tau_lag1m_CoSAR * np.conj(R_tau_lag1m)) - + dop_phase + h_phase) + + plt.figure() + csl.plot_R_tau(az, abs(R_tau_lag0), h_phase, + az, R_tau_lag0_CoSAR_a, R_tau_lag0_CoSAR_p, dB=dB, + title=r'$R_{\tau}(x,\tau=0)$', fontsize=fontsize) + if type(fileroot) == str: + filename = fileroot+'_R_tau0.eps' + plt.savefig(filename) + filename = fileroot+'_R_tau0.png' + plt.savefig(filename) + + plt.figure() + csl.plot_R_tau(az, abs(R_tau_lag1), dop_phase + h_phase, + az, R_tau_lag1_CoSAR_a, R_tau_lag1_CoSAR_p, dB=dB, + title=r'$R_{\tau}(x,\tau=1/PRF)$', fontsize=fontsize) + if type(fileroot) == str: + filename = fileroot+'_R_tau1.eps' + plt.savefig(filename) + filename = fileroot+'_R_tau1.png' + plt.savefig(filename) + + plt.figure() + csl.plot_R_tau(az, abs(R_tau_lag1m), -1.0 * dop_phase + h_phase, + az, R_tau_lag1m_CoSAR_a, R_tau_lag1m_CoSAR_p, dB=dB, + title=r'$R_{\tau}(x,\tau=-1/PRF)$', fontsize=fontsize) + if type(fileroot) == str: + filename = fileroot+'_R_tau-1.eps' + plt.savefig(filename) + filename = fileroot+'_R_tau-1.png' + plt.savefig(filename) + + +def test_geosim1d(): + """ + Test routine for geosim1d + """ + #Define geometry and radar + cs = csl.CoSARGeometry(0.8, 25, 3, 10e9, 10) + #Define ocean characteristics + sf = csl.CoSARSurface(0.1, 10, 1, 0, 1000) + #Run simulation + sim_res = geosim1d(cs, sf, 250.0, n_rg=100, PRF=10.0, Tint=300.0) + + diff --git a/cosar/geosim1d_test.py b/cosar/geosim1d_test.py new file mode 100644 index 0000000..63c93d3 --- /dev/null +++ b/cosar/geosim1d_test.py @@ -0,0 +1,50 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Jan 10 10:29:26 2014 + +@author: Paco López Dekker +""" +#%% +from cosar.geosim1d import geosim1d_analyze, geosim1d, geosim1d_montecarlo,\ + geosim1d_save, geosim1d_load +import cosar.simplesim as csl +import os + +#%% +n_it = 20 +v_stdev = 1.0 # Small to avoid big ambiguity issues +h_stdev = 1.0 +dv = 6.0 +PRF = 200.0 +tau = 0.01 +Tint = 30.0 + +respath = "/Users/plopezdekker/Documents/WORK/CoSAR/RESULTS/Figures" +os.chdir(respath) +#Define geometry and radar +cs = csl.CoSARGeometry(0.8, 25, dv, 10e9, 10) +#Define ocean characteristics +sf = csl.CoSARSurface(tau, 10, v_stdev, h_stdev, 1000) +#Run simulations +Tint = 100 +sim_res = geosim1d_montecarlo(cs, sf, 250.0, n_rg=10, PRF=PRF, Tint=Tint, + surfaceseed=1, n_it=n_it) +geosim1d_save(sim_res, 'sim_PRF200_Tint100_nr10') +sim_res = geosim1d_load('sim_PRF200_Tint100_nr10') +geosim1d_analyze(cs, sf, sim_res, fileroot='sim_PRF200_Tint100_nr10', dB=True) + +Tint = 300 +#sim_res = geosim1d_montecarlo(cs, sf, 250.0, n_rg=10, PRF=PRF, Tint=Tint, +# surfaceseed=1, n_it=n_it) +#geosim1d_save(sim_res, 'sim_PRF200_Tint300_nr10') +sim_res = geosim1d_load('sim_PRF200_Tint300_nr10') +geosim1d_analyze(cs, sf, sim_res, fileroot='sim_PRF200_Tint300_nr10', dB=True) + +Tint = 600 +#sim_res = geosim1d_montecarlo(cs, sf, 250.0, n_rg=10, PRF=PRF, Tint=Tint, +# surfaceseed=1, n_it=n_it) +#geosim1d_save(sim_res, 'sim_PRF200_Tint600_nr10') +sim_res = geosim1d_load('sim_PRF200_Tint600_nr10') +geosim1d_analyze(cs, sf, sim_res, fileroot='sim_PRF200_Tint600_nr10', dB=True) + +# %% diff --git a/cosar/results.py b/cosar/results.py new file mode 100644 index 0000000..0c1f32b --- /dev/null +++ b/cosar/results.py @@ -0,0 +1,54 @@ + # -*- coding: utf-8 -*- + +import sys +import os +import time +import numpy as np +import cosar as cosar +#from scipy import weave +from scipy import signal +import matplotlib.pyplot as plt + + +def igarss2013(figs,respath=0): + """Generate figures for IGARSS""" + if type(respath) != str: + respath = "/Users/plopezdekker/Documents/WORK/CoSAR/RESULTS/Figures" + os.chdir(respath) + if figs == 1: + filename = "sar_coh_pt_64rg_T40.png" + cosar.simplesim(256,64,tau=-1,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename) + filename = "sar_tau20_pt_64rg_T40.png" + cosar.simplesim(256,64,tau=20,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename) + filename = "sar_tau10_pt_64rg_T40.png" + cosar.simplesim(256,64,tau=10,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename) + filename = "sar_tau5_pt_64rg_T40.png" + cosar.simplesim(256,64,tau=5,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename) + filename = "sar_tau1_pt_64rg_T40.png" + cosar.simplesim(256,64,tau=1,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename) + filename = "sar_tau01_pt_64rg_T40.png" + cosar.simplesim(256,64,tau=0.1,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename) + if figs == 2: + na = 256 + nr = 256 + Tint = 40.0 + dt = 0 + taus = [20,10,5,1,0.1] + filenames = ["sar_vs_cosar_tau20_pt_256rg_T40.png", + "sar_vs_cosar_tau10_pt_256rg_T40.png", + "sar_vs_cosar_tau5_pt_256rg_T40.png", + "sar_vs_cosar_tau1_pt_256rg_T40.png", + "sar_vs_cosar_tau01_pt_256rg_T40.png"] + for ind in range(0,len(taus)): + cosar.simplesim(na,nr,tau=taus[ind],Tint=Tint,dt=dt,show_SAR=True,show_COSAR=True,dop_rms=0,pt=True,filesave=filenames[ind]) + if figs == 3: + filename = "sar_vs_cosar_tau01_seed1_1024rg_T20.png" + cosar.simplesim(256,1024,tau=0.01,Tint=20.,dt=0,show_SAR=True,show_COSAR=True,dop_rms=0.5,seed=1,filesave=filename) + nrs = [128,256,512,1024] + filenames = ["cosar_tau01_seed1_128rg_T20.png", + "cosar_tau01_seed1_256rg_T20.png", + "cosar_tau01_seed1_512rg_T20.png", + "cosar_tau01_seed1_1024rg_T20.png"] + for ind in range(0,len(nrs)): + cosar.simplesim(256,nrs[ind],tau=0.01,Tint=20.,dt=0,show_SAR=False,show_COSAR=True,dop_rms=0.5,seed=1,filesave=filenames[ind]) + \ No newline at end of file diff --git a/cosar/simplesim.py b/cosar/simplesim.py new file mode 100644 index 0000000..30e84dd --- /dev/null +++ b/cosar/simplesim.py @@ -0,0 +1,445 @@ +''' +Created on May 31, 2013 + +@author: lope_fr +''' +import sys +import time +import numpy as np +#from scipy import weave +from scipy import signal +import matplotlib.pyplot as plt +from scipy import constants +from drama import utils +from drama import geo as sar_geo +#import IPython +import matplotlib +matplotlib.rcParams['ps.fonttype'] = 42 +#matplotlib.rcParams['text.usetex'] = True + +class CoSARGeometry: + """A class containing the relevant CoSAR orbital and geometricparameters""" + _orbit_height = 35786e3 + + def __init__(self, dkr, theta_i, dv, f0, antenna_length): + self.dkr = 1.0*dkr + self.theta_i = theta_i*np.pi/180 + _dv = 1.0*np.array(dv) + if _dv.size == 1: + _dv = np.array([dv, 0.0, 0.0]) + self.dv = _dv + self.f0 = 1.0*f0 + self.L_antenna = 1.0 * antenna_length + + def wavelength(self): + return constants.c/self.f0 + + def dtheta_i(self): + res = self.dkr * self.wavelength() / (4 * np.pi * np.cos(self.theta_i)) + return res + + def h_amb(self): + #res = (2 * np.pi * np.sin(self.theta_i) * + # np.cos(self.theta_i) / self.dkr) + res = (2 * np.pi * np.cos(self.theta_i) / + (np.sin(self.theta_i) * self.dkr)) + return res + + def kz(self): + return 2 * np.pi / self.h_amb() + + def R0(self): + return sar_geo.inc_to_sr(self.theta_i, self._orbit_height) + + def L_illuminated(self): + return self.wavelength()/self.L_antenna * self.R0() + + def dDoppler_bandwidth(self): + return 2*self.dv[0]/self.L_antenna + + +class CoSARSurface: + """A class defining the scene""" + def __init__(self, tau, NRCS_stdev, vr_stdev, h_stdev, + L_NRCS, L_vr=0, L_h=0, + NRCS_mean=-10.0, vr_mean=0.0): + self.tau = tau + self.NRCS_stdev = NRCS_stdev + self.NRCS_mean = NRCS_mean + self.vr_stdev = vr_stdev + self.vr_mean = vr_stdev + self.h_stdev = h_stdev + self.L_NRCS = L_NRCS + if L_vr == 0: + self.L_vr = self.L_NRCS + else: + self.L_vr = L_vr + if L_h == 0: + self.L_h = 5.0 * self.L_NRCS + else: + self.L_h = L_h + + def scene(self, L, dx, f0): + """Generate a realization of the scene""" + naz = int(np.round(L/dx)) + az = np.linspace(-naz/2*dx, naz/2*dx, naz) + wavelength = constants.c/f0 + dop_stdev = 2*self.vr_stdev/wavelength + dop_mean = 2*self.vr_mean/wavelength + nrcs = get_nrcs(naz, int(np.round(self.L_NRCS/dx)), pt=False, + NRCS_mean=self.NRCS_mean, NRCS_stdev=self.NRCS_stdev) + dop = get_Doppler(naz, int(np.round(self.L_vr/dx)), dop_stdev) + dop = dop + dop_mean + #Reuse get_Doppler function to generate height profile + h = get_Doppler(naz, int(np.round(self.L_h/dx)), self.h_stdev) + nrcs.reshape(naz, 1) + dop.reshape(naz, 1) + h.reshape(naz, 1) + return (nrcs, dop, h, az) + + def L_min(self): + """Returns the minimum length scale""" + return np.min(np.array([self.L_h, self.L_NRCS, self.L_vr])) + + +def Q_w(R, NRCS, dx, N_r, T_i, tau_ca, az_res, F_n=1.0): + """ + Estimate CoSAR imaging quality. + R is the space-dependent expected correlation + NRCS is the corresponding NRCS + dx is the grid spacing + N_r the number of range gates averaged + T_i the integration time + tau_ca the decorrelation time + az_res the azimuth resolution + """ + Q_scaling = N_r / F_n * T_i / tau_ca * az_res**2 + Power = np.sum(NRCS) * dx + Q_w = Q_scaling * np.abs(R)**2 / Power**2 + return Q_w + + +def get_nrcs(naz, gaussian_length, pt=False, SCR=30, + NRCS_mean=-10, NRCS_stdev=10.0): + """Generate a random but smooth NRCS""" + #We use a uniform distribution + dist_amp = np.sqrt(12) * NRCS_stdev + nrcs_db = dist_amp * np.random.random(naz) - dist_amp / 2 + NRCS_mean + #A Gaussian Filter + gfilt = signal.gaussian(int(naz / 4), gaussian_length, sym=True) + gfilt = gfilt / gfilt.sum() + nrcs_db_filt = signal.fftconvolve(nrcs_db, gfilt, mode='same') + nrcs_filt = utils.db2lin(nrcs_db_filt) + if pt: + nrcs_filt[:] = utils.db2lin(-SCR) + nrcs_filt[[naz / 4, naz / 2]] = 1 + return nrcs_filt + + +def get_Doppler(naz, gaussian_length, dop_rms=1): + """Generate random Doppler frequencies""" + dop = np.random.randn(naz) + #A Gaussian Filter + gfilt = signal.gaussian(int(naz / 4), gaussian_length, sym=True) + gfilt = gfilt / gfilt.sum() + dop_filt = signal.fftconvolve(dop, gfilt, mode='same') + norm_factor = dop_rms/np.std(dop_filt) + dop_filt = dop_filt * norm_factor + return dop_filt + + +def cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0, han=False, dt=0): + """CoSAR processing assuming no range migrations""" + #Cross correlate signals + if dt == 0: + corr_est = np.mean(rx2 * np.conjugate(rx1), axis=1) + else: + if dt > 0: + corr_p2p = rx2[dt:, :] * np.conjugate(rx1)[:-dt, :] + else: + corr_p2p = rx2[0:dt, :] * np.conjugate(rx1)[-dt:, :] + corr_est = np.mean(corr_p2p, axis=1) + + naz = az.size + nt = corr_est.size + if han: + az_win = np.hanning(nt) + else: + az_win = np.ones(nt) + + T_i_scaling = np.mean(az_win**2) + #Result + s_est = np.zeros(naz, dtype=np.complex) + for t_ind in range(0, corr_est.size): + #Calculate ranges from targets to both radars + #r_2 = ((pos1[t_ind, 1])**2 + (az - pos1[t_ind, 0])**2 + + # pos1[t_ind, 2] ** 2 ) + r1 = np.sqrt((pos1[t_ind, 1])**2 + (az - pos1[t_ind, 0])**2 + + pos1[t_ind, 2]**2) + + r2 = np.sqrt((pos2[t_ind + dt, 1])**2 + + (az - pos2[t_ind + dt, 0])**2 + + pos2[t_ind+dt, 2]**2) + + dphase = 2*k0*(r2-r1) + phasor = np.exp(1j*dphase) + s_est = s_est + az_win[t_ind] * phasor * corr_est[t_ind] + #Normalize + s_est = s_est # / corr_est.size / 2 + return (s_est, T_i_scaling) + + +def sar_proc_1d(rx, pos, az, k0, han=False): + """SAR processing assuming no range migrations using backprojection""" + nt = rx.shape[0] + nrg = rx.shape[1] + naz = az.size + if han: + az_win = np.hanning(nt) + else: + az_win = np.ones(nt) + #Results + T_i_scaling = np.mean(az_win**2) + s_est = np.zeros([naz, nrg], dtype=np.complex) + for t_ind in range(0, nt): + r = (np.sqrt((pos[t_ind, 1])**2 + (az - pos[t_ind, 0])**2 + + pos[t_ind, 2]**2)) + phase = 2 * k0 * r + phasor = np.exp(1j * phase) + s_est = s_est + az_win[t_ind] * rx[t_ind, :] * phasor.reshape([naz, 1]) + pow_est = np.mean(np.abs(s_est)**2, axis=1) + return (pow_est/nt, T_i_scaling) + + +def plot_nrcs(az_true, nrcs_true, az_proc, nrcs_proc, dB=True, nrcs_err=0): + """Routine to plot cosar estimated NRCS and ground truth""" + ax = plt.subplot(1, 1, 1) + res_dim = nrcs_proc.shape + if np.size(res_dim) == 2: + n_res = res_dim[0] + else: + n_res = 1 + nrcs_proc.shape = [n_res, res_dim[0]] + if dB: + p_true_nrcs, = plt.plot(az_true, utils.db(nrcs_true, norm=True), 'b', + label="True") + p_cosar_nrcs, = plt.plot(az_proc, utils.db(nrcs_proc[0, :], norm=True), + 'r', label="CoSAR") + for it_res in range(1, n_res): + plt.plot(az_proc, utils.db(nrcs_proc[it_res, :], norm=True), 'r') + if np.size(nrcs_err) == np.size(nrcs_true): + p_error, = plt.plot(az_true, + utils.db(nrcs_true + nrcs_err, norm=True), + 'b--') + + plt.xlabel('Azimuth [km]') + plt.ylabel('NRCS [dB]') + else: + p_true_nrcs, = plt.plot(az_true, nrcs_true, 'b', + label="True") + p_cosar_nrcs, = plt.plot(az_proc, nrcs_proc[0, :], 'r', + label="CoSAR") + for it_res in range(1, n_res): + plt.plot(az_proc, nrcs_proc[it_res, :], 'r') + if np.size(nrcs_err) == np.size(nrcs_true): + p_error, = plt.plot(az_true, + nrcs_true + nrcs_err, 'b--') + plt.xlabel('Azimuth [km]') + plt.ylabel('NRCS') + #Add labels + handles, labels = ax.get_legend_handles_labels() + plt.legend(handles, labels, loc=3) + plt.show() + + +def plot_R_tau(az_true, r_true_a, r_true_p, + az_proc, r_proc_a, r_proc_p, dB=True, r_err=0, + title=False, fontsize=12): + """ + Routine to plot CoSAR estimated R_tau and ground truth + Variables ending in _a are amplitude, those ending in _p are phases + """ + ax = plt.subplot(2, 1, 1) + res_dim = r_proc_a.shape + if np.size(res_dim) == 2: + n_res = res_dim[0] + else: + n_res = 1 + r_proc_a.shape = [n_res, res_dim[0]] + r_proc_p.shape = [n_res, res_dim[0]] + if dB: + p_true_r, = plt.plot(az_true, utils.db(r_true_a, norm=True), 'b', + label="True") + p_cosar_r, = plt.plot(az_proc, + utils.db(r_proc_a[0, :], norm=True), + 'r', label="CoSAR") + for it_res in range(1, n_res): + plt.plot(az_proc, + utils.db(r_proc_a[it_res, :], norm=True), + 'r') + if np.size(r_err) == np.size(r_true_a): + p_error, = plt.plot(az_true, + utils.db(r_true_a + r_err, norm=True), + 'b--') + + #plt.xlabel('Azimuth [km]') + plt.ylabel(r'$|R|$ $\mathrm{[dB]}$', fontsize=fontsize) + else: + p_true_r, = plt.plot(az_true, r_true_a, 'b', + label="True") + p_cosar_r, = plt.plot(az_proc, r_proc_a[0, :], 'r', + label="CoSAR") + for it_res in range(1, n_res): + plt.plot(az_proc, r_proc_a[it_res, :], 'r') + if np.size(r_err) == np.size(r_true_a): + p_error, = plt.plot(az_true, + r_true_a + r_err, 'b--') + #plt.xlabel('Azimuth [km]') + plt.ylabel(r'$|R|$', fontsize=fontsize) + if type(title) == str: + plt.title(title, fontsize=fontsize) + #Add labels + handles, labels = ax.get_legend_handles_labels() + plt.legend(handles, labels, loc=3, fontsize=fontsize) + #Now plot the phase + ax = plt.subplot(2, 1, 2) + p_true_r_ph, = plt.plot(az_true, r_true_p, 'b', label='True') + p_true_r_ph, = plt.plot(az_proc, r_proc_p[0, :], 'r', label='CoSAR') + for it_res in range(1, n_res): + plt.plot(az_proc, r_proc_p[it_res, :], 'r') + + plt.xlabel('Azimuth [km]', fontsize=fontsize) + plt.ylabel(r'$\mathrm{arg}(R)$ $\mathrm{[rad]}$', fontsize=fontsize) + + + +def simplesim_show(az,true_nrcs,est_nrcs,true_Doppler=0,est_Doppler=0,sar_nrcs=0, + show_SAR=False,show_Doppler=False,show_COSAR=True,filesave=0): + """Routine to plot simplesim simulation results""" + if show_Doppler: + ax = plt.subplot(2,1,1) + else: + ax = plt.subplot(1,1,1) + p_true_nrcs, = plt.plot(az,utils.db(true_nrcs,norm=True),'b',label="True") + if show_COSAR: + p_cosar_nrcs, = plt.plot(az,utils.db(est_nrcs,norm=True),'r',label="CoSAR") + plt.xlabel('Azimuth') + plt.ylabel('Intensity [dB]') + if show_SAR: + plt.plot(az,utils.db(sar_nrcs,norm=True),'k--',label="SAR") + #Add labels + handles, labels = ax.get_legend_handles_labels() + plt.legend(handles,labels,loc=3) + if show_Doppler > 0: + plt.subplot(2,1,2) + plt.plot(az,true_Doppler,'b',az,est_Doppler,'r') + plt.xlabel('Azimuth') + plt.ylabel('Doppler [Hz]') + plt.show() + #Check if we want to print this + if type(filesave) == str: + plt.savefig(filesave) + plt.close() + +def simplesim_scene_init(naz,length_scale,dop_rms=1.0,pt='False'): + nrcs = get_nrcs(naz,length_scale,pt=pt) + dop = get_Doppler(naz,10,dop_rms) + nrcs.reshape(naz,1) + dop.reshape(naz,1) + return (nrcs,dop) + +def simplesim(naz,nrg,tau=0.0,pt=False,hamming=True,dop_rms=1.0,Tint=20.0,dt=0, + show_SAR=False,show_COSAR=True,filesave=0,seed=None): + + #Create some complexarray of random numbers + + dx = 1.0 + #Spacecraft height and ground range + z1 = z2 = 5e3 + x1 = x2 = -5e3 + R0 = 10e3 + v1, v2 = 1.0, -1.0 + f0 = 10e9 + k0 = 2*np.pi*f0/constants.c + l0 = constants.c/f0 + fDop_max = naz * dx / R0 * 2 * abs(v2-v1)/l0 + PRF = 8*fDop_max + print("PRF = %i", int(PRF)) + t_v = np.linspace(-Tint/2, Tint/2, Tint*PRF) + az = np.linspace(-naz/2*dx, naz/2*dx, naz) + ##Seed simulator if needed + np.random.seed(seed) + #nrcs + nrcs, dop = simplesim_scene_init(naz, 10, pt=pt) + dop_phasor = np.exp(2j*np.pi*dop/PRF) + #print dop_phasor.shape + expected_amp = np.sqrt(nrcs) + #teporal decorrelation factor or whatever + # Get trajectory of spacecraft + pos1 = np.zeros([Tint*PRF, 3]) + pos2 = np.zeros([Tint*PRF, 3]) + pos1[:, 0] = v1*t_v + pos2[:, 0] = v2*t_v + pos1[:, 1] = x1 + pos2[:, 1] = x2 + pos1[:, 2] = z1 + pos2[:, 2] = z2 + #Correlation coefficcient + # alpha**(tau*PRF) = 1/e + # tau*PRF * log(alpha)=-1 alpha = exp(-1/(tau*PRF) + alpha = 1 + if tau > 0: + alpha = np.exp(-1.0/(tau*PRF)) + if nrg > 1: + uscat = np.random.randn(naz, nrg) + 1j*np.random.randn(naz, nrg) + else: + uscat = np.exp(1j*2*np.pi*np.random.rand(naz, nrg)) + for t_ind in range(0, t_v.size): + #Generate random complex scatterers + uscat = (alpha * uscat * dop_phasor.reshape(naz, 1) + + np.sqrt(1-alpha**2) * (np.random.randn(naz, nrg) + + 1j * np.random.randn(naz, nrg))) + #Multipliy them with a fixed amplitude + scat = uscat * expected_amp.reshape(naz, 1) + #Calculate ranges from targets to both radars + r1 = ((np.sqrt((pos1[t_ind, 1])**2 + + (az - pos1[t_ind, 0])**2 + pos1[t_ind, 2]**2)).reshape(naz, 1)) + r2 = ((np.sqrt((pos2[t_ind, 1])**2 + + (az - pos2[t_ind, 0])**2 + pos2[t_ind, 2]**2)).reshape(naz, 1)) + #Include 2-way phase in field + s1 = scat * np.exp(-2j*k0*r1) + s2 = scat * np.exp(-2j*k0*r2) + # Compute range profiles by integrating in azimuth. + # This assumed that there is no range migration + rx1_ = s1.sum(axis=0) + rx2_ = s2.sum(axis=0) + # Save in a temporal stack + if t_ind == 0: + rx1 = rx1_.copy() + rx2 = rx2_.copy() + else: + rx1 = np.vstack((rx1, rx1_)) + rx2 = np.vstack((rx2, rx2_)) + #Now we have a range-slow time stack of CoSAR data + proc_img, T_i_scaling = cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0, + han=hamming) + if dt > 0: + proc_ximg = cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0, han=hamming, + dt=dt) + #The phase of this should be the Doppler phase + est_Doppler = np.angle(proc_ximg)/(dt/PRF)/(2*np.pi) + else: + est_Doppler = 0 + + if show_SAR: + sar_img = sar_proc_1d(rx1, pos1, az, k0, han=hamming) + else: + sar_img = 0 + #IPython.embed() + #Show results + #plt.ion() + simplesim_show(az, nrcs, np.abs(proc_img), dop, est_Doppler, sar_img, + show_SAR=show_SAR, show_Doppler=(dt > 0), + show_COSAR=show_COSAR, filesave=filesave) + print("Done!") + diff --git a/cosar/viciprep.py b/cosar/viciprep.py new file mode 100644 index 0000000..e646ae7 --- /dev/null +++ b/cosar/viciprep.py @@ -0,0 +1,67 @@ +__author__ = "Paco Lopez Dekker" +__email__ = "F.LopezDekker@tudeft.nl" + +import numpy as np +import scipy.constants as cnst + +def spectral_shift(f0, theta_i, dtheta_i): + return f0 * dtheta_i / np.tan(theta_i) + + +def CoSAR_airborne(target_res, h=2.5e3, Bhor=50, v_sar = 100, dv_sar=2, theta_i_deg=45, az_bmwdth_deg=30, + f0=np.array([1.24e9, 5.4e9, 9.65e9]), + Bw=np.array([200, 200, 500]) * 1e6, + antenna_length=None): + + + wl0 = cnst.c / f0 + theta_i = np.radians(theta_i_deg) + U = 5 + tau_c = 3 * wl0 / U + Bn = np.cos(theta_i) * Bhor + R0 = h / np.cos(theta_i) + az_bmwdth = np.radians(az_bmwdth_deg) + if not antenna_length is None: + az_bmwdth = wl0 / antenna_length + + # Spectral shift + dtheta_i = Bn / R0 + df = spectral_shift(f0, theta_i, dtheta_i) + print("Spectral shift [MHz]") + print(df/1e6) + #az_bmwdth = np.radians(25) + + + width_gr = R0 * np.sin(az_bmwdth) + + Lap = wl0 * R0 / 2 / target_res + T_CoSAR_i = Lap / dv_sar + + ## CoSAR Quality factor + Bw_common = Bw - df + dgr = cnst.c / 2 / Bw_common / np.sin(theta_i) + Nr = target_res / dgr + + Q = Nr * T_CoSAR_i / tau_c * target_res ** 2 / az_bmwdth + + # Printing of the resuls + print("Slant range: %f" % R0) + + np.set_printoptions(precision=3) + print("Azimuth antenna foot print:") + print(width_gr) + print("SAR focused res (decorrelation limited)") + print(wl0 / (v_sar * tau_c)* R0/2) + print("CoSAR aperture lengths") + print(Lap) + print("Ground Distance travelled") + print(T_CoSAR_i * v_sar) + print("CoSAR integation time") + print(T_CoSAR_i) + print("Number of range gates") + print(Nr.astype(np.int)) + print("Quality factor (Q):") + print(Q) + print("Pseudo SCR (dB):") + print(5*np.log10(Q)) + np.set_printoptions() \ No newline at end of file diff --git a/cosarsim.py b/cosarsim.py new file mode 100644 index 0000000..6e46e91 --- /dev/null +++ b/cosarsim.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python + +''' +COSAR Simulator launcher + +Usage :: + ./cosar (GUI version) + ./cosar -nogui /path/to/parfile (Console version) + +''' + +import os +import sys +import argparse +import cosar + +# Check Python version +if tuple(sys.version_info[:2]) < (2, 7): + raise Exception('Python %d.%d is not supported. Please use version 2.6 ' + 'or greater.\n' % tuple(sys.version_info[:2])) + +def main(argv): + """Launches COSAR Simulator""" + #Handle command line argument + parser = argparse.ArgumentParser(description="Simple CoSAR simulation") + parser.add_argument("-p", "--point_target", action="store_true") + parser.add_argument("-w", "--hamming", action="store_true") + parser.add_argument("-a", "--azimuth_points", type=int, default = 256) + parser.add_argument("-r", "--range_samples", type=int, default = 1024) + parser.add_argument("-c", "--correlation_coefficient", type=float,default = 0.0) + args = parser.parse_args() + naz = args.azimuth_points + nrg = args.range_samples + alpha = args.correlation_coefficient + cosar.simplesim(naz, nrg, alpha=alpha, pt=args.point_target, hamming=args.hamming) + + + +if __name__ == '__main__': + main(sys.argv) \ No newline at end of file diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..3fc7a49 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,27 @@ +[metadata] +name = CoSAR +version = 0.1 +author = Paco López-Dekker +author_email = F.LopezDekker@tudelft.nl +description = CoSAR simulation code +long_description = file: README.rst +long_description_content_type = text/x-rst +# url = https://gitlab.tudelft.nl/plopezdekker/s1sea +#project_urls = +# Bug Tracker = https://gitlab.tudelft.nl/plopezdekker/waddensar/-/issues +classifiers = + Programming Language :: Python :: 3 + License :: OSI Approved :: GNU General Public License v3 (GPLv3) + Operating System :: OS Independent + Intended Audience :: Science/Research + Topic :: Scientific/Engineering + +[options] +packages = find: +python_requires = >=3.9 +install_requires = + numpy>=1.20 + scipy + numexpr>=2.7 + matplotlib + drama diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..d50db6d --- /dev/null +++ b/setup.py @@ -0,0 +1,15 @@ +from setuptools import setup, find_packages + +setup( + name='cosar', + version='24.03.26', + # packages=['oceansar'], + packages=find_packages(exclude=['contrib', 'docs', 'tests*']), + # package_dir={'': 'oceansar'}, + # url='url = https://gitlab.tudelft.nl/plopezdekker/s1sea', + license='GPL-3.0', + author='Paco Lopez Dekker', + author_email='F.LopezDekker@tudelft.nl', + description='', + install_requires=['numpy', 'scipy', 'matplotlib', 'drama'] +) -- GitLab