diff --git a/src/mintpy/cli/view.py b/src/mintpy/cli/view.py
index caaf44969..1df824629 100755
--- a/src/mintpy/cli/view.py
+++ b/src/mintpy/cli/view.py
@@ -20,6 +20,7 @@
   view.py velocity.h5 --ref-yx  210 566                              #change reference pixel for display
   view.py velocity.h5 --sub-lat 31.05 31.10 --sub-lon 130.05 130.10  #subset in lalo / yx
   view.py velocity.h5 velocity --mask waterBody.h5 --mask-vmax 1
+  view.py velocity.h5 velocity --style scatter --scatter-size 12
 
   view.py timeseries.h5
   view.py timeseries.h5 --ref-date 20101120     #change reference date
@@ -79,6 +80,12 @@ def create_parser(subparsers=None):
                              '  reverse = x * -1\n'
                              '  inverse = 1 / x')
 
+    # plot data in different styles: image, scatter, contour etc.
+    parser.add_argument('--style', dest='style', choices={'image', 'scatter'}, default='image',
+                        help='Plot data as image or scatter (default: %(default)s).')
+    parser.add_argument('--scatter-size', dest='scatter_marker_size', type=float, metavar='SIZE', default=10,
+                        help='Scatter marker size in points**2 (default: %(default)s).')
+
     parser = arg_utils.add_data_disp_argument(parser)
     parser = arg_utils.add_dem_argument(parser)
     parser = arg_utils.add_figure_argument(parser)
diff --git a/src/mintpy/view.py b/src/mintpy/view.py
index aee29a3b8..bc4fd7fe2 100644
--- a/src/mintpy/view.py
+++ b/src/mintpy/view.py
@@ -476,6 +476,18 @@ def plot_slice(ax, data, metadata, inps):
     global vprint
     vprint = print if inps.print_msg else lambda *args, **kwargs: None
 
+    def extent2meshgrid(extent: tuple, ds_shape: list):
+        """Get mesh grid coordinates for a given extent and shape.
+        Parameters: extent - tuple of float for (left, right, bottom, top) in data coordinates
+                    shape  - list of int for [length, width] of the data
+        Returns:    xx/yy  - 1D np.ndarray of the data coordinates
+        """
+        height, width = ds_shape
+        x = np.linspace(extent[0], extent[1], width)
+        y = np.linspace(extent[2], extent[3], height)[::-1]  # reverse the Y-axis
+        xx, yy = np.meshgrid(x, y)
+        return xx.flatten(), yy.flatten()
+
     # colormap: str -> object
     if isinstance(inps.colormap, str):
         inps.colormap = pp.ColormapExt(
@@ -497,6 +509,9 @@ def plot_slice(ax, data, metadata, inps):
     num_row, num_col = data.shape
     lalo_digit = ut.get_lalo_digit4display(metadata, coord_unit=inps.coord_unit)
 
+    # common options for data visualization
+    kwargs = dict(cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1], alpha=inps.transparency, zorder=1)
+
     #----------------------- Plot in Geo-coordinate --------------------------------------------#
     if (inps.geo_box
             and inps.coord_unit.startswith(('deg', 'meter'))
@@ -539,11 +554,20 @@ def plot_slice(ax, data, metadata, inps):
         # Plot data
         if inps.disp_dem_blend:
             im = pp.plot_blend_image(ax, data, dem, inps, print_msg=inps.print_msg)
+
+        elif inps.style == 'image':
+            vprint(f'plotting data as {inps.style} via matplotlib.pyplot.imshow ...')
+            im = ax.imshow(data, extent=inps.extent, origin='upper', interpolation=inps.interpolation,
+                           animated=inps.animation, **kwargs)
+
+        elif inps.style == 'scatter':
+            vprint(f'plotting data as {inps.style} via matplotlib.pyplot.scatter (can take some time) ...')
+            xx, yy = extent2meshgrid(inps.extent, data.shape)
+            im = ax.scatter(xx, yy, c=data.flatten(), marker='o', s=inps.scatter_marker_size, **kwargs)
+            ax.axis('equal')
+
         else:
-            vprint('plotting data ...')
-            im = ax.imshow(data, cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1],
-                           extent=inps.extent, origin='upper', interpolation=inps.interpolation,
-                           alpha=inps.transparency, animated=inps.animation, zorder=1)
+            raise ValueError(f'Un-recognized plotting style: {inps.style}!')
 
         # Draw faultline using GMT lonlat file
         if inps.faultline_file:
@@ -657,11 +681,19 @@ def format_coord(x, y):
         # Plot Data
         if inps.disp_dem_blend:
             im = pp.plot_blend_image(ax, data, dem, inps, print_msg=inps.print_msg)
+
+        elif inps.style == 'image':
+            vprint('plotting data via matplotlib.pyplot.imshow ...')
+            im = ax.imshow(data, extent=inps.extent, interpolation=inps.interpolation, **kwargs)
+
+        elif inps.style == 'scatter':
+            vprint('plotting data via matplotlib.pyplot.scatter (can take some time) ...')
+            xx, yy = extent2meshgrid(inps.extent, data.shape)
+            im = ax.scatter(xx, yy, c=data.flatten(), marker='o', s=inps.scatter_marker_size, **kwargs)
+            ax.axis('equal')
+
         else:
-            vprint('plotting data ...')
-            im = ax.imshow(data, cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1],
-                           extent=inps.extent, interpolation=inps.interpolation,
-                           alpha=inps.transparency, zorder=1)
+            raise ValueError(f'Un-recognized plotting style: {inps.style}!')
         ax.tick_params(labelsize=inps.font_size)
 
         # Plot Seed Point
@@ -783,11 +815,11 @@ def format_coord(x, y):
     # rotate Y-axis tick labels
     # link: https://stackoverflow.com/questions/10998621
     if inps.ylabel_rot:
-        kwargs = dict(rotation=inps.ylabel_rot)
+        tick_kwargs = dict(rotation=inps.ylabel_rot)
         # center the vertical alignment for vertical tick labels
         if inps.ylabel_rot % 90 == 0:
-            kwargs['va'] = 'center'
-        plt.setp(ax.get_yticklabels(), **kwargs)
+            tick_kwargs['va'] = 'center'
+        plt.setp(ax.get_yticklabels(), **tick_kwargs)
         vprint(f'rotate Y-axis tick labels by {inps.ylabel_rot} deg')
 
     return ax, inps, im, cbar