Dirty beam images added!

app.py

@@ -29,6 +29,8 @@ import dash_core_components as dcc
 import dash_html_components as html
 import dash_bootstrap_components as dbc
 import plotly.express as px
+from plotly.subplots import make_subplots
+import plotly.graph_objects as go
 from astropy.time import Time
 from astropy import coordinates as coord
 from astropy import units as u
@@ -827,11 +829,10 @@ def main_page(results_visible=False, summary_output=None, fig_elev_output=None,
                         html.Br(),
                         html.Div([
                             html.Br(),
-                            html.H4("Resulting dirty map"),
+                            html.H4("Resulting dirty beam"),
                             html.Br()
                         ]),
-                        html.Div(children=[dcc.Graph(id='fig-dirtymap', figure=fig_dirty_map_output,
-                                            relayoutData={'xaxis': {'autorange': 'reverse'}})],
+                        html.Div(children=[dcc.Graph(id='fig-dirtymap', figure=fig_dirty_map_output)],
                                  className='tex2jax_ignore')
                     ])])
                 ]),
@@ -1355,10 +1356,30 @@ def get_fig_uvplane(obs):
 
 
 def get_fig_dirty_map(obs):
-    dirty_map, laxis = obs.get_dirtymap(pixsize=1024)
-    print(laxis)
-    return px.imshow(img=dirty_map, x=laxis, y=laxis, labels={'x': 'RA (mas)', 'y': 'Dec (mas)'}, \
+    dirty_map_nat, laxis = obs.get_dirtymap(pixsize=1024, robust='natural')
+    dirty_map_uni, laxis = obs.get_dirtymap(pixsize=1024, robust='uniform')
+    fig1 = px.imshow(img=dirty_map_nat, x=laxis, y=laxis[::-1], labels={'x': 'RA (mas)', 'y': 'Dec (mas)'}, \
             aspect='equal')
+    fig2 = px.imshow(img=dirty_map_uni, x=laxis, y=laxis[::-1], labels={'x': 'RA (mas)', 'y': 'Dec (mas)'}, \
+            aspect='equal')
+    # for i,f in enumerate([fig1, fig2]):
+    #     f.layout.xaxis.autorange = "reversed"
+    #     f.layout.yaxis.autorange = True
+    #     f['layout']['coloraxis']['showscale'] = False
+
+    # fig1.layout.xaxis.autorange = "reversed"
+    # fig2.layout.xaxis.autorange = "reversed"
+    fig = make_subplots(rows=1, cols=2, subplot_titles=('Natural weighting', 'Uniform weighting'),
+                        shared_xaxes=True, shared_yaxes=True)
+    fig.add_trace(fig1.data[0], row=1, col=1)
+    fig.add_trace(fig2.data[0], row=1, col=2)
+    fig.update_layout(coloraxis={'showscale': False, 'colorscale': 'Inferno'}, showlegend=False, xaxis={'autorange': "reversed"},
+                      yaxis={'autorange': True}, xaxis2={'autorange': "reversed"}, autosize=False)
+    fig.update_xaxes(title_text="RA (mas)", constrain="domain")
+    fig.update_yaxes(title_text="Dec (mas)", row=1, col=1, scaleanchor="x", scaleratio=1)
+
+    return fig
+    return fig1
             # layout :  xaxis autorange reversed
 
 

vlbiplanobs/graphical_elements.py

@@ -210,7 +210,7 @@ def summary_card_beam(app, obs):
     temp_msg += [html.Br(), html.P([f"The expected synthesized beam will be approx. {synthbeam['bmaj'].to(synthbeam_units).value:.3n} x {synthbeam['bmin'].to(synthbeam_units):.3n}", html.Sup("2"), \
             f", PA = {synthbeam['pa']:.3n}."])]
     temp_msg += [html.P("Note that the synthesized beam can significantly change depending "
-                        "on the weighting used during imaging (null weighting assumed here).")]
+                        "on the weighting used during imaging (natural weighting assumed here).")]
     return create_sensitivity_card('Resolution', temp_msg)
 
 
@@ -321,7 +321,7 @@ def summary_card_rms(app, obs):
                 html.Div(className='row justify-content-center',
                 children=html.Img(src=app.get_asset_url("waves.png"), width='100%',
                                   height='75rem', style={'display': 'inline'}))]
-    temp_msg += [html.P(f"The expected rms thermal noise for your target is {rms:.3n}/beam when no weighting is applied during imaging. Note that ~20% higher values may be expected for RFI-contaminated bands.")]
+    temp_msg += [html.P(f"The expected rms thermal noise for your target is {rms:.3n}/beam using natural weighting during imaging. Note that ~20% higher values may be expected for RFI-contaminated bands.")]
     temp_msg += [html.P(f"The achieved sensitivity implies a rms of {rms_channel:.3n}/beam per spectral "
                         f"channel, or approx. {rms_time:.3n}/beam per time integration "
                         f"({optimal_units(obs.inttime, [u.s,u.ms,u.us]):.3n}).")]

vlbiplanobs/observation.py

@@ -841,11 +841,15 @@ class Observation(object):
         return self._synth_beam
 
 
-    def get_dirtymap(self, pixsize:int =1024):
+    def get_dirtymap(self, pixsize:int = 1024, robust: str = "natural"):
         """Returns the dirty beam produced for the given observation.
 
-        Input:  - pixsize : int
+        Input:
+            - pixsize : int
                 Size in pixels of the returned (squared) image. By default 1024.
+            - robust : str
+                The weighting Briggs robust used to compute the dirty map.
+                It must be either 'natural' or 'uniform'.
 
         Returns:
             - dirty_image : (pixsize x pixsize) np.array
@@ -854,18 +858,23 @@ class Observation(object):
                 An array representing the values of each pixel in the image in mas for each axis.
 
         """
+        assert robust in ('natural', 'uniform')
         pixsize = 1024
         uvimg = np.zeros((pixsize, pixsize))
 
         # Gridding uv.  Inspired from VNSIM code (https://github.com/ZhenZHAO/VNSIM)
         uvscaling = pixsize/(0.9*self.longest_baseline()[1].to(u.cm).value)
         uvdata = self.get_uv_array()
-        for uv in uvdata:
-            uvimg[int(pixsize//2 + round(uv[0]*uvscaling)), int(pixsize//2 +round(uv[1]*uvscaling))] = 1
+        if robust == 'natural':
+            for uv in uvdata:
+                uvimg[int(pixsize//2 + round(uv[0]*uvscaling)), int(pixsize//2 +round(uv[1]*uvscaling))] += 1
+        else:
+            for uv in uvdata:
+                uvimg[int(pixsize//2 + round(uv[0]*uvscaling)), int(pixsize//2 +round(uv[1]*uvscaling))] = 1
 
-        dirty_beam = np.real(np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(uvimg))))
+        dirty_beam = np.real(np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(uvimg/np.max(uvimg)))))
         imgsize = (uvscaling*u.rad/2).to(u.mas) # angular equivalent size of the resulting image
-        return dirty_beam/np.max(dirty_beam), np.linspace(-imgsize, imgsize, pixsize)
+        return dirty_beam.T/np.max(dirty_beam), np.linspace(-imgsize, imgsize, pixsize)