diff --git a/applications/topsApp.py b/applications/topsApp.py
index 6dd2b65..66baa8b 100755
--- a/applications/topsApp.py
+++ b/applications/topsApp.py
@@ -519,6 +519,22 @@ ION_AZSHIFT_FLAG = Application.Parameter('ION_azshiftFlag',
     mandatory=False,
     doc='')
 
+#seperated islands or areas usually affect ionosphere estimation and it's better to mask them
+#out. check ion/ion_cal/raw_no_projection.ion for areas to be masked out.
+#The parameter is a 2-D list. Each element in the 2-D list is a four-element list: [firstLine,
+#lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the
+#four elements is specified as -1, the program will use firstLine/lastLine/firstColumn/
+#lastColumn instead. For exmple, if you want to mask the following two areas out, you can
+#specify a 2-D list like:
+#[[100, 200, 100, 200],[1000, 1200, 500, 600]]
+ION_MASKED_AREAS = Application.Parameter('ION_maskedAreas',
+    public_name = 'areas masked out in ionospheric phase estimation',
+    default = None,
+    type = int,
+    mandatory = False,
+    container = list,
+    doc = 'areas masked out in ionospheric phase estimation')
+
 ION_NUMBER_AZIMUTH_LOOKS = Application.Parameter('ION_numberAzimuthLooks',
     public_name='total number of azimuth looks in the ionosphere processing',
     default=50,
@@ -674,6 +690,7 @@ class TopsInSAR(Application):
                       ION_IONSHIFT_FILTERING_WINSIZE_MAX,
                       ION_IONSHIFT_FILTERING_WINSIZE_MIN,
                       ION_AZSHIFT_FLAG,
+                      ION_MASKED_AREAS,
                       ION_NUMBER_AZIMUTH_LOOKS,
                       ION_NUMBER_RANGE_LOOKS,
                       ION_NUMBER_AZIMUTH_LOOKS0,
diff --git a/components/isceobj/TopsProc/runIon.py b/components/isceobj/TopsProc/runIon.py
index 06b48a1..c81655f 100644
--- a/components/isceobj/TopsProc/runIon.py
+++ b/components/isceobj/TopsProc/runIon.py
@@ -14,6 +14,7 @@ import isceobj
 import logging
 from isceobj.Constants import SPEED_OF_LIGHT
 from isceobj.TopsProc.runBurstIfg import loadVirtualArray
+from isceobj.Alos2Proc.runIonFilt import reformatMaskedAreas
 
 
 logger = logging.getLogger('isce.topsinsar.ion')
@@ -89,6 +90,7 @@ def setup(self):
     ionParam.ionshiftFilteringWinsizeMax = self.ION_ionshiftFilteringWinsizeMax
     ionParam.ionshiftFilteringWinsizeMin = self.ION_ionshiftFilteringWinsizeMin
     ionParam.azshiftFlag = self.ION_azshiftFlag
+    ionParam.maskedAreas = self.ION_maskedAreas
 
     ionParam.numberAzimuthLooks = self.ION_numberAzimuthLooks
     ionParam.numberRangeLooks = self.ION_numberRangeLooks
@@ -1283,6 +1285,14 @@ def ionosphere(self, ionParam):
     cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
     amp = np.sqrt(lowerAmp**2+upperAmp**2)
 
+    #masked out user-specified areas
+    if ionParam.maskedAreas != None:
+        maskedAreas = reformatMaskedAreas(ionParam.maskedAreas, length, width)
+        for area in maskedAreas:
+            lowerUnw[area[0]:area[1], area[2]:area[3]] = 0
+            upperUnw[area[0]:area[1], area[2]:area[3]] = 0
+            cor[area[0]:area[1], area[2]:area[3]] = 0
+
     #compute ionosphere
     fl = SPEED_OF_LIGHT / ionParam.radarWavelengthLower
     fu = SPEED_OF_LIGHT / ionParam.radarWavelengthUpper
@@ -1549,6 +1559,14 @@ def ionSwathBySwath(self, ionParam):
         cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
         amp = np.sqrt(lowerAmp**2+upperAmp**2)
 
+        #masked out user-specified areas
+        if ionParam.maskedAreas != None:
+            maskedAreas = reformatMaskedAreas(ionParam.maskedAreas, length, width)
+            for area in maskedAreas:
+                lowerUnw[area[0]:area[1], area[2]:area[3]] = 0
+                upperUnw[area[0]:area[1], area[2]:area[3]] = 0
+                cor[area[0]:area[1], area[2]:area[3]] = 0
+
         #compute ionosphere
         fl = SPEED_OF_LIGHT / ionParam.radarWavelengthLower
         fu = SPEED_OF_LIGHT / ionParam.radarWavelengthUpper
@@ -1923,6 +1941,13 @@ def filt_gaussian(self, ionParam):
     cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
     amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
 
+    #masked out user-specified areas
+    if ionParam.maskedAreas != None:
+        maskedAreas = reformatMaskedAreas(ionParam.maskedAreas, length, width)
+        for area in maskedAreas:
+            ion[area[0]:area[1], area[2]:area[3]] = 0
+            cor[area[0]:area[1], area[2]:area[3]] = 0
+
     ########################################################################################
     #AFTER COHERENCE IS RESAMPLED AT grd2ion, THERE ARE SOME WIRED VALUES
     cor[np.nonzero(cor<0)] = 0.0
diff --git a/examples/input_files/topsApp.xml b/examples/input_files/topsApp.xml
index ae3241d..cc96c81 100644
--- a/examples/input_files/topsApp.xml
+++ b/examples/input_files/topsApp.xml
@@ -200,6 +200,16 @@ By default, ESD is always performed.
 #2: use full burst
 <property name="correct phase error caused by ionosphere azimuth shift">1</property>
 
+#seperated islands or areas usually affect ionosphere estimation and it's better to mask them
+#out. check ion/ion_cal/raw_no_projection.ion for areas to be masked out.
+#The parameter is a 2-D list. Each element in the 2-D list is a four-element list: [firstLine,
+#lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the
+#four elements is specified as -1, the program will use firstLine/lastLine/firstColumn/
+#lastColumn instead. For exmple, if you want to mask the following two areas out, you can
+#specify a 2-D list like:
+#[[100, 200, 100, 200],[1000, 1200, 500, 600]]
+<property name="areas masked out in ionospheric phase estimation">None</property>
+
 #better NOT try changing the following two parameters, since they are related
 #to the filtering parameters above
 <property name="total number of azimuth looks in the ionosphere processing">50</property>