diff --git a/pyproject.toml b/pyproject.toml
index fd2eeff50e87247643d943f275d6daac00f475d9..b476a6f3b2b70582be5b0b39c6819052167f91b2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -45,6 +45,7 @@ dependencies = [
     "python-casacore",
     "python-dateutil",
     "reproject",
+    "requests",
     "Rtree",
     "scipy",
     "shapely",
diff --git a/rapthor/lib/facet.py b/rapthor/lib/facet.py
index 065c8a0ea2d6dca1ea16d2d86f262b14a19d10cd..b627e95b50b4e83219bc0415ed03de1fb367281b 100644
--- a/rapthor/lib/facet.py
+++ b/rapthor/lib/facet.py
@@ -2,8 +2,225 @@
 Module that holds functions and classes related to faceting
 """
 import numpy as np
-import scipy as sp
 import scipy.spatial
+from shapely.geometry import Point, Polygon
+from shapely.prepared import prep
+from PIL import Image, ImageDraw
+from astropy.coordinates import Angle
+import ast
+import logging
+from rapthor.lib import miscellaneous as misc
+from matplotlib import patches
+import lsmtool
+from lsmtool import tableio
+import tempfile
+
+
+class Facet(object):
+    """
+    Base Facet class
+
+    Parameters
+    ----------
+    name : str
+        Name of facet
+    ra : float or str
+        RA of reference coordinate in degrees (if float) or as a string in a
+        format supported by astropy.coordinates.Angle
+    dec : float or str
+        Dec of reference coordinate in degrees (if float) or as a string in a
+        format supported by astropy.coordinates.Angle
+    vertices : list of tuples
+        List of (RA, Dec) tuples, one for each vertex of the facet
+    """
+    def __init__(self, name, ra, dec, vertices):
+        self.name = name
+        self.log = logging.getLogger('rapthor:{0}'.format(self.name))
+        if type(ra) is str:
+            ra = Angle(ra).to('deg').value
+        if type(dec) is str:
+            dec = Angle(dec).to('deg').value
+        self.ra = misc.normalize_ra(ra)
+        self.dec = misc.normalize_dec(dec)
+        self.vertices = np.array(vertices)
+
+        # Convert input (RA, Dec) vertices to (x, y) polygon
+        self.wcs = make_wcs(self.ra, self.dec)
+        self.polygon_ras = [radec[0] for radec in self.vertices]
+        self.polygon_decs = [radec[1] for radec in self.vertices]
+        x_values, y_values = radec2xy(self.wcs, self.polygon_ras, self.polygon_decs)
+        polygon_vertices = [(x, y) for x, y in zip(x_values, y_values)]
+        self.polygon = Polygon(polygon_vertices)
+
+        # Find the size and center coordinates of the facet
+        xmin, ymin, xmax, ymax = self.polygon.bounds
+        self.size = min(0.5, max(xmax-xmin, ymax-ymin) *
+                        abs(self.wcs.wcs.cdelt[0]))  # degrees
+        self.x_center = xmin + (xmax - xmin)/2
+        self.y_center = ymin + (ymax - ymin)/2
+        self.ra_center, self.dec_center = xy2radec(self.wcs, self.x_center, self.y_center)
+
+    def set_skymodel(self, skymodel):
+        """
+        Sets the facet's sky model
+
+        The input sky model is filtered to contain only those sources that lie
+        inside the facet's polygon. The filtered sky model is stored in
+        self.skymodel
+
+        Parameters
+        ----------
+        skymodel : LSMTool skymodel object
+            Input sky model
+        """
+        self.skymodel = filter_skymodel(self.polygon, skymodel, self.wcs)
+
+    def download_panstarrs(self, max_search_cone_radius=0.5):
+        """
+        Returns a Pan-STARRS sky model for the area around the facet
+
+        Note: the resulting sky model may contain sources outside the facet's
+        polygon
+
+        Parameters
+        ----------
+        max_search_cone_radius : float, optional
+            The maximum radius in degrees to use in the cone search. The smaller
+            of this radius and the minimum radius that covers the facet is used
+
+        Returns
+        -------
+        skymodel : LSMTool skymodel object
+            The Pan-STARRS sky model
+        """
+        try:
+            with tempfile.NamedTemporaryFile() as fp:
+                misc.download_skymodel(self.ra_center, self.dec_center, fp.name,
+                                       radius=min(max_search_cone_radius, self.size/2),
+                                       overwrite=True, source='PANSTARRS')
+                skymodel = lsmtool.load(fp.name)
+                skymodel.group('every')
+        except IOError:
+            # Comparison catalog not downloaded successfully
+            self.log.warning('The Pan-STARRS catalog could not be successfully '
+                             'downloaded')
+            skymodel = tableio.makeEmptyTable()
+
+        return skymodel
+
+    def find_astrometry_offsets(self, comparison_skymodel=None, min_number=5):
+        """
+        Finds the astrometry offsets for sources in the facet
+
+        The offsets are calculated as (LOFAR model value) - (comparison model
+        value); e.g., a positive Dec offset indicates that the LOFAR sources
+        are on average North of the comparison source positions.
+
+        The offsets are stored in self.astrometry_diagnostics, a dict with
+        the following keys (see LSMTool's compare operation for details of the
+        diagnostics):
+
+            'meanRAOffsetDeg', 'stdRAOffsetDeg', 'meanClippedRAOffsetDeg',
+            'stdClippedRAOffsetDeg', 'meanDecOffsetDeg', 'stdDecOffsetDeg',
+            'meanClippedDecOffsetDeg', 'stdClippedDecOffsetDeg'
+
+        Note: if the comparison is unsuccessful, self.astrometry_diagnostics is
+        an empty dict
+
+        Parameters
+        ----------
+        comparison_skymodel : LSMTool skymodel object, optional
+            Comparison sky model. If not given, the Pan-STARRS catalog is
+            used
+        min_number : int, optional
+            Minimum number of sources required for comparison
+        """
+        self.astrometry_diagnostics = {}
+        if comparison_skymodel is None:
+            comparison_skymodel = self.download_panstarrs()
+
+        # Find the astrometry offsets between the facet's sky model and the
+        # comparison sky model
+        #
+        # Note: If there are no successful matches, the compare() method
+        # returns None
+        if len(comparison_skymodel) >= min_number:
+            result = self.skymodel.compare(comparison_skymodel,
+                                           radius='5 arcsec',
+                                           excludeMultiple=True,
+                                           make_plots=False)
+            # Save offsets
+            if result is not None:
+                self.astrometry_diagnostics.update(result)
+        else:
+            self.log.warning('Too few matches to determine astrometry offsets '
+                             '(min_number = {0} but number of matches '
+                             '= {1})'.format(min_number, len(comparison_skymodel)))
+
+    def get_matplotlib_patch(self, wcs=None):
+        """
+        Returns a matplotlib patch for the facet polygon
+
+        Parameters
+        ----------
+        wcs : WCS object, optional
+            WCS object defining (RA, Dec) <-> (x, y) transformation. If not given,
+            the facet's transformation is used
+
+        Returns
+        -------
+        patch : matplotlib patch object
+            The patch for the facet polygon
+        """
+        if wcs is not None:
+            x, y = radec2xy(wcs, self.polygon_ras, self.polygon_decs)
+        else:
+            x = self.polygon.exterior.coords.xy[0]
+            y = self.polygon.exterior.coords.xy[1]
+        xy = np.vstack([x, y]).transpose()
+        patch = patches.Polygon(xy=xy, edgecolor='black', facecolor='white')
+
+        return patch
+
+
+class SquareFacet(Facet):
+    """
+    Wrapper class for a square facet
+
+    Parameters
+    ----------
+    name : str
+        Name of facet
+    ra : float or str
+        RA of reference coordinate in degrees (if float) or as a string in a
+        format supported by astropy.coordinates.Angle
+    dec : float or str
+        Dec of reference coordinate in degrees (if float) or as a string in a
+        format supported by astropy.coordinates.Angle
+    width : float
+        Width in degrees of facet
+    """
+    def __init__(self, name, ra, dec, width):
+        if type(ra) is str:
+            ra = Angle(ra).to('deg').value
+        if type(dec) is str:
+            dec = Angle(dec).to('deg').value
+        ra = misc.normalize_ra(ra)
+        dec = misc.normalize_dec(dec)
+        wcs = make_wcs(ra, dec)
+
+        # Make the vertices
+        xmin = wcs.wcs.crpix[0] - width / 2 / abs(wcs.wcs.cdelt[0])
+        xmax = wcs.wcs.crpix[0] + width / 2 / abs(wcs.wcs.cdelt[0])
+        ymin = wcs.wcs.crpix[1] - width / 2 / abs(wcs.wcs.cdelt[1])
+        ymax = wcs.wcs.crpix[1] + width / 2 / abs(wcs.wcs.cdelt[1])
+        ra_llc, dec_llc = xy2radec(wcs, xmin, ymin)  # (RA, Dec) of lower-left corner
+        ra_tlc, dec_tlc = xy2radec(wcs, xmin, ymax)  # (RA, Dec) of top-left corner
+        ra_trc, dec_trc = xy2radec(wcs, xmax, ymax)  # (RA, Dec) of top-right corner
+        ra_lrc, dec_lrc = xy2radec(wcs, xmax, ymin)  # (RA, Dec) of lower-right corner
+        vertices = [(ra_llc, dec_llc), (ra_tlc, dec_tlc), (ra_trc, dec_trc), (ra_lrc, dec_lrc)]
+
+        super().__init__(name, ra, dec, vertices)
 
 
 def make_facet_polygons(ra_cal, dec_cal, ra_mid, dec_mid, width_ra, width_dec):
@@ -38,13 +255,13 @@ def make_facet_polygons(ra_cal, dec_cal, ra_mid, dec_mid, width_ra, width_dec):
     if width_ra <= 0.0 or width_dec <= 0.0:
         raise ValueError('The RA/Dec width cannot be zero or less')
     wcs_pixel_scale = 20.0 / 3600.0  # 20"/pixel
-    wcs = makeWCS(ra_mid, dec_mid, wcs_pixel_scale)
+    wcs = make_wcs(ra_mid, dec_mid, wcs_pixel_scale)
     x_cal, y_cal = radec2xy(wcs, ra_cal, dec_cal)
-    x_mid, y_mid = radec2xy(wcs, [ra_mid], [dec_mid])
+    x_mid, y_mid = radec2xy(wcs, ra_mid, dec_mid)
     width_x = width_ra / wcs_pixel_scale / 2.0
     width_y = width_dec / wcs_pixel_scale / 2.0
-    bounding_box = np.array([x_mid[0] - width_x, x_mid[0] + width_x,
-                             y_mid[0] - width_y, y_mid[0] + width_y])
+    bounding_box = np.array([x_mid - width_x, x_mid + width_x,
+                             y_mid - width_y, y_mid + width_y])
 
     # Tessellate and convert resulting facet polygons from (x, y) to (RA, Dec)
     vor = voronoi(np.stack((x_cal, y_cal)).T, bounding_box)
@@ -56,13 +273,13 @@ def make_facet_polygons(ra_cal, dec_cal, ra_mid, dec_mid, width_ra, width_dec):
         facet_polys.append(vertices)
     facet_points = []
     for point in vor.filtered_points:
-        ra, dec = xy2radec(wcs, [point[0]], [point[1]])
-        facet_points.append((ra[0], dec[0]))
+        ra, dec = xy2radec(wcs, point[0], point[1])
+        facet_points.append((ra, dec))
 
     return facet_points, facet_polys
 
 
-def radec2xy(wcs, RA, Dec):
+def radec2xy(wcs, ra, dec):
     """
     Returns x, y for input RA, Dec
 
@@ -70,57 +287,97 @@ def radec2xy(wcs, RA, Dec):
     ----------
     wcs : WCS object
         WCS object defining transformation
-    RA : list
-        List of RA values in degrees
-    Dec : list
-        List of Dec values in degrees
+    ra : float, list, or numpy array
+        RA value(s) in degrees
+    dec : float, list, or numpy array
+        Dec value(s) in degrees
 
     Returns
     -------
-    x, y : list, list
-        Lists of x and y pixel values corresponding to the input RA and Dec
+    x, y : float, list, or numpy array
+        x and y pixel values corresponding to the input RA and Dec
         values
     """
-    x = []
-    y = []
-
-    for ra_deg, dec_deg in zip(RA, Dec):
+    x_list = []
+    y_list = []
+    if type(ra) is list or type(ra) is np.ndarray:
+        ra_list = ra
+    else:
+        ra_list = [float(ra)]
+    if type(dec) is list or type(dec) is np.ndarray:
+        dec_list = dec
+    else:
+        dec_list = [float(dec)]
+    if len(ra_list) != len(dec_list):
+        raise ValueError('RA and Dec must be of equal length')
+
+    for ra_deg, dec_deg in zip(ra_list, dec_list):
         ra_dec = np.array([[ra_deg, dec_deg]])
-        x.append(wcs.wcs_world2pix(ra_dec, 0)[0][0])
-        y.append(wcs.wcs_world2pix(ra_dec, 0)[0][1])
+        x_list.append(wcs.wcs_world2pix(ra_dec, 0)[0][0])
+        y_list.append(wcs.wcs_world2pix(ra_dec, 0)[0][1])
+
+    # Return the same type as the input
+    if type(ra) is list or type(ra) is np.ndarray:
+        x = x_list
+    else:
+        x = x_list[0]
+    if type(dec) is list or type(dec) is np.ndarray:
+        y = y_list
+    else:
+        y = y_list[0]
     return x, y
 
 
 def xy2radec(wcs, x, y):
     """
-    Returns input RA, Dec for input x, y
+    Returns RA, Dec for input x, y
 
     Parameters
     ----------
     wcs : WCS object
         WCS object defining transformation
-    x : list
-        List of x values in pixels
-    y : list
-        List of y values in pixels
+    x : float, list, or numpy array
+        x value(s) in pixels
+    y : float, list, or numpy array
+        y value(s) in pixels
 
     Returns
     -------
-    RA, Dec : list, list
-        Lists of RA and Dec values corresponding to the input x and y pixel
+    RA, Dec : float, list, or numpy array
+        RA and Dec values corresponding to the input x and y pixel
         values
     """
-    RA = []
-    Dec = []
-
-    for xp, yp in zip(x, y):
+    ra_list = []
+    dec_list = []
+    if type(x) is list or type(x) is np.ndarray:
+        x_list = x
+    else:
+        x_list = [float(x)]
+    if type(y) is list or type(y) is np.ndarray:
+        y_list = y
+    else:
+        y_list = [float(y)]
+    if len(x_list) != len(y_list):
+        raise ValueError('x and y must be of equal length')
+
+    for xp, yp in zip(x_list, y_list):
         x_y = np.array([[xp, yp]])
-        RA.append(wcs.wcs_pix2world(x_y, 0)[0][0])
-        Dec.append(wcs.wcs_pix2world(x_y, 0)[0][1])
-    return RA, Dec
-
-
-def makeWCS(ra, dec, wcs_pixel_scale=10.0/3600.0):
+        ra_list.append(wcs.wcs_pix2world(x_y, 0)[0][0])
+        dec_list.append(wcs.wcs_pix2world(x_y, 0)[0][1])
+
+    # Return the same type as the input
+    if type(x) is list or type(x) is np.ndarray:
+        ra = ra_list
+    else:
+        ra = ra_list[0]
+    if type(y) is float or type(y) is np.ndarray:
+        dec = dec_list
+    else:
+        dec = dec_list[0]
+    return ra, dec
+
+
+def make_wcs(ra, dec, wcs_pixel_scale=10.0/3600.0):
     """
     Makes simple WCS object
 
@@ -214,7 +471,7 @@ def voronoi(cal_coords, bounding_box):
 
     # Compute Voronoi, sorting the output regions to match the order of the
     # input coordinates
-    vor = sp.spatial.Voronoi(points)
+    vor = scipy.spatial.Voronoi(points)
     sorted_regions = np.array(vor.regions, dtype=object)[np.array(vor.point_region)]
     vor.regions = sorted_regions.tolist()
 
@@ -229,8 +486,8 @@ def voronoi(cal_coords, bounding_box):
             else:
                 x = vor.vertices[index, 0]
                 y = vor.vertices[index, 1]
-                if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
-                       bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
+                if not (bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
+                        bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
                     flag = False
                     break
         if region and flag:
@@ -280,3 +537,120 @@ def make_ds9_region_file(center_coords, facet_polygons, outfile, names=None):
 
     with open(outfile, 'w') as f:
         f.writelines(lines)
+
+
+def read_ds9_region_file(region_file):
+    """
+    Read a ds9 facet region file and return facets
+
+    Parameters
+    ----------
+    region_file : str
+        Filename of input ds9 region file
+
+    Returns
+    -------
+    facets : list
+        List of Facet objects
+    """
+    facets = []
+    with open(region_file, 'r') as f:
+        lines = f.readlines()
+    for line in lines:
+        # Each facet in the region file is defined by two consecutive lines:
+        #   - the first starts with 'polygon' and gives the (RA, Dec) vertices
+        #   - the second starts with 'point' and gives the reference (RA, Dec)
+        #     and the facet name
+        if line.startswith('polygon'):
+            vertices = ast.literal_eval(line.split('polygon')[1])
+            polygon_ras = [ra for ra in vertices[::2]]
+            polygon_decs = [dec for dec in vertices[1::2]]
+            vertices = [(ra, dec) for ra, dec in zip(polygon_ras, polygon_decs)]
+        if line.startswith('point'):
+            ra, dec = ast.literal_eval(line.split('point')[1])
+            if 'text' in line:
+                name = line.split('text=')[1].strip()
+            else:
+                name = f'facet_{ra}_{dec}'
+            facets.append(Facet(name, ra, dec, vertices))
+
+    return facets
+
+
+def filter_skymodel(polygon, skymodel, wcs):
+    """
+    Filters input skymodel to select only sources that lie inside the input facet
+
+    Parameters
+    ----------
+    polygon : Shapely polygon object
+        Polygon object to use for filtering
+    skymodel : LSMTool skymodel object
+        Input sky model to be filtered
+    wcs : WCS object
+        WCS object defining image to sky transformations
+
+    Returns
+    -------
+    filtered_skymodel : LSMTool skymodel object
+        Filtered sky model
+    """
+    # Make list of sources
+    RA = skymodel.getColValues('Ra')
+    Dec = skymodel.getColValues('Dec')
+    x, y = radec2xy(wcs, RA, Dec)
+    x = np.array(x)
+    y = np.array(y)
+
+    # Keep only those sources inside the bounding box
+    inside = np.zeros(len(skymodel), dtype=bool)
+    xmin, ymin, xmax, ymax = polygon.bounds
+    inside_ind = np.where((x >= xmin) & (x <= xmax) & (y >= ymin) & (y <= ymax))
+    inside[inside_ind] = True
+    skymodel.select(inside)
+    if len(skymodel) == 0:
+        return skymodel
+    RA = skymodel.getColValues('Ra')
+    Dec = skymodel.getColValues('Dec')
+    x, y = radec2xy(wcs, RA, Dec)
+    x = np.array(x)
+    y = np.array(y)
+
+    # Now check the actual boundary against filtered sky model. We first do a quick (but
+    # coarse) check using ImageDraw with a padding of at least a few pixels to ensure the
+    # quick check does not remove sources spuriously. We then do a slow (but precise)
+    # check using Shapely
+    xpadding = max(int(0.1 * (max(x) - min(x))), 3)
+    ypadding = max(int(0.1 * (max(y) - min(y))), 3)
+    xshift = int(min(x)) - xpadding
+    yshift = int(min(y)) - ypadding
+    xsize = int(np.ceil(max(x) - min(x))) + 2*xpadding
+    ysize = int(np.ceil(max(y) - min(y))) + 2*ypadding
+    x -= xshift
+    y -= yshift
+    prepared_polygon = prep(polygon)
+
+    # Unmask everything outside of the polygon + its border (outline)
+    inside = np.zeros(len(skymodel), dtype=bool)
+    mask = Image.new('L', (xsize, ysize), 0)
+    verts = [(xv-xshift, yv-yshift) for xv, yv in zip(polygon.exterior.coords.xy[0],
+                                                      polygon.exterior.coords.xy[1])]
+    ImageDraw.Draw(mask).polygon(verts, outline=1, fill=1)
+    inside_ind = np.where(np.array(mask).transpose()[(x.astype(int), y.astype(int))])
+    inside[inside_ind] = True
+
+    # Now check sources in the border precisely
+    mask = Image.new('L', (xsize, ysize), 0)
+    ImageDraw.Draw(mask).polygon(verts, outline=1, fill=0)
+    border_ind = np.where(np.array(mask).transpose()[(x.astype(int), y.astype(int))])
+    points = [Point(xs, ys) for xs, ys in zip(x[border_ind], y[border_ind])]
+    indexes = []
+    for i in range(len(points)):
+        indexes.append(border_ind[0][i])
+    i_points = zip(indexes, points)
+    i_outside_points = [(i, p) for (i, p) in i_points if not prepared_polygon.contains(p)]
+    for idx, _ in i_outside_points:
+        inside[idx] = False
+    skymodel.select(inside)
+
+    return skymodel
diff --git a/rapthor/lib/miscellaneous.py b/rapthor/lib/miscellaneous.py
index a832379658c4f83c1a02e9211d2febc1c6caef67..2e05f31b3fabbf4f6ded7c258d18f58dd102b923 100644
--- a/rapthor/lib/miscellaneous.py
+++ b/rapthor/lib/miscellaneous.py
@@ -21,38 +21,34 @@ import lsmtool
 from scipy.interpolate import interp1d
 import astropy.units as u
 import mocpy
+import requests
 
 
 def download_skymodel(ra, dec, skymodel_path, radius=5.0, overwrite=False, source='TGSS',
                       targetname='Patch'):
     """
-    Download the skymodel for the target field
+    Downloads a skymodel for the given position and radius
 
     Parameters
     ----------
     ra : float
-        Right ascension of the skymodel centre.
+        Right ascension in degrees of the skymodel centre
     dec : float
-        Declination of the skymodel centre.
+        Declination in degrees of the skymodel centre
     skymodel_path : str
-        Full name (with path) to the skymodel.
+        Full name (with path) to the output skymodel
     radius : float, optional
-        Radius for the TGSS/GSM cone search in degrees.
+        Radius for the cone search in degrees. For Pan-STARRS, the radius must be
+        <= 0.5 degrees
     source : str, optional
-        Source where to obtain a skymodel from. Can be TGSS or GSM. Default is TGSS.
+        Source where to obtain a skymodel from. Can be one of: TGSS, GSM, LOTSS, or
+        PANSTARRS. Note: the PANSTARRS sky model is only suitable for use in
+        astrometry checks and should not be used for calibration
     overwrite : bool, optional
-        Overwrite the existing skymodel pointed to by skymodel_path.
+        Overwrite the existing skymodel pointed to by skymodel_path
     target_name : str, optional
-        Give the patch a certain name. Default is "Patch".
+        Give the patch a certain name
     """
-    SKY_SERVERS = {'TGSS': 'http://tgssadr.strw.leidenuniv.nl/cgi-bin/gsmv4.cgi?'
-                           'coord={ra:f},{dec:f}&radius={radius:f}&unit=deg&deconv=y',
-                   'GSM': 'https://lcs165.lofar.eu/cgi-bin/gsmv1.cgi?'
-                          'coord={ra:f},{dec:f}&radius={radius:f}&unit=deg&deconv=y',
-                   'LOTSS': ''}  # Server is empty since we handle this through LSMTool.
-    if source.upper() not in SKY_SERVERS.keys():
-        raise ValueError('Unsupported sky model source specified! Please use LOTSS, TGSS or GSM.')
-
     logger = logging.getLogger('rapthor:skymodel')
 
     file_exists = os.path.isfile(skymodel_path)
@@ -77,22 +73,31 @@ def download_skymodel(ra, dec, skymodel_path, radius=5.0, overwrite=False, sourc
                        'existing sky model!'.format(skymodel_path))
         os.remove(skymodel_path)
 
-    # Check if LoTSS has coverage.
-    if source.upper().strip() == 'LOTSS':
+    # Check the radius for Pan-STARRS (it must be <= 0.5 degrees)
+    source = source.upper().strip()
+    if source == 'PANSTARRS' and radius > 0.5:
+        raise ValueError('The radius for Pan-STARRS must be <= 0.5 deg')
+
+    # Check if LoTSS has coverage
+    if source == 'LOTSS':
         logger.info('Checking LoTSS coverage for the requested centre and radius.')
         mocpath = os.path.join(os.path.dirname(skymodel_path), 'dr2-moc.moc')
-        subprocess.run(['wget', 'https://lofar-surveys.org/public/DR2/catalogues/dr2-moc.moc', '-O', mocpath], capture_output=True, check=True)
+        subprocess.run(['wget', 'https://lofar-surveys.org/public/DR2/catalogues/dr2-moc.moc',
+                        '-O', mocpath], capture_output=True, check=True)
         moc = mocpy.MOC.from_fits(mocpath)
         covers_centre = moc.contains(ra * u.deg, dec * u.deg)
-        # Checking single coordinates, so get rid of the array.
+
+        # Checking single coordinates, so get rid of the array
         covers_left = moc.contains(ra * u.deg - radius * u.deg, dec * u.deg)[0]
         covers_right = moc.contains(ra * u.deg + radius * u.deg, dec * u.deg)[0]
         covers_bottom = moc.contains(ra * u.deg, dec * u.deg - radius * u.deg)[0]
         covers_top = moc.contains(ra * u.deg, dec * u.deg + radius * u.deg)[0]
         if covers_centre and not (covers_left and covers_right and covers_bottom and covers_top):
-            logger.warning('Incomplete LoTSS coverage for the requested centre and radius! Please check the field coverage in plots/field_coverage.png!')
+            logger.warning('Incomplete LoTSS coverage for the requested centre and radius! '
+                           'Please check the field coverage in plots/field_coverage.png!')
         elif not covers_centre and (covers_left or covers_right or covers_bottom or covers_top):
-            logger.warning('Incomplete LoTSS coverage for the requested centre and radius! Please check the field coverage in plots/field_coverage.png!')
+            logger.warning('Incomplete LoTSS coverage for the requested centre and radius! '
+                           'Please check the field coverage in plots/field_coverage.png!')
         elif not covers_centre and not (covers_left and covers_right and covers_bottom and covers_top):
             raise ValueError('No LoTSS coverage for the requested centre and radius!')
         else:
@@ -101,32 +106,55 @@ def download_skymodel(ra, dec, skymodel_path, radius=5.0, overwrite=False, sourc
     logger.info('Downloading skymodel for the target into ' + skymodel_path)
     max_tries = 5
     for tries in range(1, 1 + max_tries):
-        if source.upper().strip() == 'LOTSS':
+        retry = False
+        if source == 'LOTSS' or source == 'TGSS' or source == 'GSM':
             try:
-                lotssmodel = lsmtool.skymodel.SkyModel('lotss', VOPosition=[ra, dec], VORadius=radius)
-                lotssmodel.write(skymodel_path)
-                if len(lotssmodel) > 0:
+                skymodel = lsmtool.skymodel.SkyModel(source, VOPosition=[ra, dec], VORadius=radius)
+                skymodel.write(skymodel_path)
+                if len(skymodel) > 0:
                     break
             except ConnectionError:
-                if tries == max_tries:
-                    raise IOError('Download of LoTSS sky model failed after {} attempts.'.format(max_tries))
-                else:
-                    logger.error('Attempt #{0:d} to download LoTSS sky model failed. Attempting '
-                                 '{1:d} more times.'.format(tries, max_tries - tries))
-                    time.sleep(5)
-        elif (source.upper().strip() == 'TGSS') or (source.upper().strip() == 'GSM'):
-            result = subprocess.run(['wget', '-O', skymodel_path,
-                                    SKY_SERVERS[source].format(ra=ra, dec=dec, radius=radius)])
-            if result.returncode != 0:
-                if tries == max_tries:
-                    raise IOError('Download of TGSS sky model failed after {} '
-                                  'attempts.'.format(max_tries))
+                retry = True
+        elif source == 'PANSTARRS':
+            baseurl = 'https://catalogs.mast.stsci.edu/api/v0.1/panstarrs'
+            release = 'dr1'  # the release with the mean data
+            table = 'mean'  # the main catalog, with the mean data
+            cat_format = 'csv'  # use csv format for the intermediate file
+            url = f'{baseurl}/{release}/{table}.{cat_format}'
+            search_params = {'ra': ra,
+                             'dec': dec,
+                             'radius': radius,
+                             'nDetections.min': '5',  # require detection in at least 5 epochs
+                             'columns': ['objID', 'ramean', 'decmean']  # get only the info we need
+                             }
+            try:
+                result = requests.get(url, params=search_params, timeout=300)
+                if result.ok:
+                    # Convert the result to makesourcedb format and write to the output file
+                    lines = result.text.split('\n')[1:]  # split and remove header line
+                    out_lines = ['FORMAT = Name, Ra, Dec, Type, I, ReferenceFrequency=1e6\n']
+                    for line in lines:
+                        # Add entries for type and Stokes I flux density
+                        if line.strip():
+                            out_lines.append(line.strip() + ',POINT,0.0,\n')
+                    with open(skymodel_path, 'w') as f:
+                        f.writelines(out_lines)
+                    break
                 else:
-                    logger.error('Attempt #{0:d} to download TGSS sky model failed. Attempting '
-                                 '{1:d} more times.'.format(tries, max_tries - tries))
-                    time.sleep(5)
+                    retry = True
+            except requests.exceptions.Timeout:
+                retry = True
+        else:
+            raise ValueError('Unsupported sky model source specified! Please use LOTSS, TGSS, '
+                             'GSM, or PANSTARRS.')
+
+        if retry:
+            if tries == max_tries:
+                raise IOError('Download of {0} sky model failed after {1} attempts.'.format(source, max_tries))
             else:
-                break
+                logger.error('Attempt #{0:d} to download {1} sky model failed. Attempting '
+                             '{2:d} more times.'.format(tries, source, max_tries - tries))
+                time.sleep(5)
 
     if not os.path.isfile(skymodel_path):
         raise IOError('Sky model file "{}" does not exist after trying to download the '
diff --git a/rapthor/lib/sector.py b/rapthor/lib/sector.py
index 480085bced4e6d70d797191af667853697971bc1..a4ce51b6594e66068440821f55d8608c9633bc4c 100644
--- a/rapthor/lib/sector.py
+++ b/rapthor/lib/sector.py
@@ -4,13 +4,11 @@ Definition of the Sector class that holds parameters for an image or predict sec
 import logging
 import numpy as np
 from rapthor.lib import miscellaneous as misc
-from rapthor.lib import cluster
+from rapthor.lib import cluster, facet
 import lsmtool
 from astropy.coordinates import Angle,  SkyCoord
 import astropy.units as u
-from shapely.geometry import Point, Polygon
-from shapely.prepared import prep
-from PIL import Image, ImageDraw
+from shapely.geometry import Polygon
 import pickle
 import os
 import copy
@@ -62,6 +60,8 @@ class Sector(object):
         self.imsize = None  # set to None to force calculation in set_imaging_parameters()
         self.wsclean_image_padding = 1.2  # the WSClean default value, used in the workflows
         self.diagnostics = []  # list to hold dicts of image diagnostics
+        self.calibration_skymodel = None  # set by Field.update_skymodel()
+        self.max_nmiter = None  # set by the strategy
 
         # Make copies of the observation objects, as each sector may have its own
         # observation-specific settings
@@ -358,61 +358,7 @@ class Sector(object):
         filtered_skymodel : LSMTool skymodel object
             Filtered sky model
         """
-        # Make list of sources
-        RA = skymodel.getColValues('Ra')
-        Dec = skymodel.getColValues('Dec')
-        x, y = self.field.radec2xy(RA, Dec)
-        x = np.array(x)
-        y = np.array(y)
-
-        # Keep only those sources inside the sector bounding box
-        inside = np.zeros(len(skymodel), dtype=bool)
-        xmin, ymin, xmax, ymax = self.poly.bounds
-        inside_ind = np.where((x >= xmin) & (x <= xmax) & (y >= ymin) & (y <= ymax))
-        inside[inside_ind] = True
-        skymodel.select(inside)
-        if len(skymodel) == 0:
-            return skymodel
-        RA = skymodel.getColValues('Ra')
-        Dec = skymodel.getColValues('Dec')
-        x, y = self.field.radec2xy(RA, Dec)
-        x = np.array(x)
-        y = np.array(y)
-
-        # Now check the actual sector boundary against filtered sky model
-        xpadding = max(int(0.1 * (max(x) - min(x))), 3)
-        ypadding = max(int(0.1 * (max(y) - min(y))), 3)
-        xshift = int(min(x)) - xpadding
-        yshift = int(min(y)) - ypadding
-        xsize = int(np.ceil(max(x) - min(x))) + 2*xpadding
-        ysize = int(np.ceil(max(y) - min(y))) + 2*ypadding
-        x -= xshift
-        y -= yshift
-        prepared_polygon = prep(self.poly)
-
-        # Unmask everything outside of the polygon + its border (outline)
-        inside = np.zeros(len(skymodel), dtype=bool)
-        mask = Image.new('L', (xsize, ysize), 0)
-        verts = [(xv-xshift, yv-yshift) for xv, yv in zip(self.poly.exterior.coords.xy[0],
-                                                          self.poly.exterior.coords.xy[1])]
-        ImageDraw.Draw(mask).polygon(verts, outline=1, fill=1)
-        inside_ind = np.where(np.array(mask).transpose()[(x.astype(int), y.astype(int))])
-        inside[inside_ind] = True
-
-        # Now check sources in the border precisely
-        mask = Image.new('L', (xsize, ysize), 0)
-        ImageDraw.Draw(mask).polygon(verts, outline=1, fill=0)
-        border_ind = np.where(np.array(mask).transpose()[(x.astype(int), y.astype(int))])
-        points = [Point(xs, ys) for xs, ys in zip(x[border_ind], y[border_ind])]
-        indexes = []
-        for i in range(len(points)):
-            indexes.append(border_ind[0][i])
-        i_points = zip(indexes, points)
-        i_outside_points = [(i, p) for (i, p) in i_points if not prepared_polygon.contains(p)]
-        for idx, _ in i_outside_points:
-            inside[idx] = False
-        skymodel.select(inside)
-        return skymodel
+        return facet.filter_skymodel(self.poly, skymodel, self.field.wcs)
 
     def get_obs_parameters(self, parameter):
         """
diff --git a/rapthor/operations/image.py b/rapthor/operations/image.py
index 7a7662f46f1980378909eefc39e91cae0b0e620c..4dcfc1c918c7c8af211b705f0611cddcf06c527a 100644
--- a/rapthor/operations/image.py
+++ b/rapthor/operations/image.py
@@ -5,6 +5,7 @@ import os
 import json
 import logging
 import shutil
+import glob
 from rapthor.lib import miscellaneous as misc
 from rapthor.lib.operation import Operation
 from rapthor.lib.cwl import CWLFile, CWLDir
@@ -307,6 +308,16 @@ class Image(Operation):
                         shutil.rmtree(dst_filename)
                     shutil.copytree(src_filename, dst_filename)
 
+            # The astrometry and photometry plots
+            dst_dir = os.path.join(self.parset['dir_working'], 'plots', 'image_{}'.format(self.index))
+            misc.create_directory(dst_dir)
+            diagnostic_plots = glob.glob(os.path.join(self.pipeline_working_dir, '*.pdf'))
+            for src_filename in diagnostic_plots:
+                dst_filename = os.path.join(dst_dir, os.path.basename(src_filename))
+                if os.path.exists(dst_filename):
+                    os.remove(dst_filename)
+                shutil.copy(src_filename, dst_filename)
+
             # Read in the image diagnostics and log a summary of them
             diagnostics_file = image_root + '.image_diagnostics.json'
             with open(diagnostics_file, 'r') as f:
diff --git a/rapthor/pipeline/parsets/image_pipeline.cwl b/rapthor/pipeline/parsets/image_pipeline.cwl
index 4a4f10573f58bf72bfac4568a863ab244ae6b900..ba83b370c61c3076b6eef3e08e952c6138ba6980 100644
--- a/rapthor/pipeline/parsets/image_pipeline.cwl
+++ b/rapthor/pipeline/parsets/image_pipeline.cwl
@@ -433,6 +433,18 @@ outputs:
     outputSource:
       - image_sector/sector_diagnostics
     type: File[]
+  - id: sector_offsets
+    outputSource:
+      - image_sector/sector_offsets
+    type: File[]
+  - id: sector_diagnostic_plots
+    outputSource:
+      - image_sector/sector_diagnostic_plots
+    type:
+      type: array
+      items:
+        type: array
+        items: File
   - id: visibilities
     outputSource:
       - image_sector/visibilities
@@ -662,6 +674,8 @@ steps:
       - id: sector_skymodels
 {% endif %}
       - id: sector_diagnostics
+      - id: sector_offsets
+      - id: sector_diagnostic_plots
 {% if use_facets %}
       - id: sector_region_file
 {% endif %}
diff --git a/rapthor/pipeline/parsets/image_sector_pipeline.cwl b/rapthor/pipeline/parsets/image_sector_pipeline.cwl
index 31185495210fd731544b2ae64a6c43b61772513e..b60270153e25f8dad90c7be76eb620f668bd4350 100644
--- a/rapthor/pipeline/parsets/image_sector_pipeline.cwl
+++ b/rapthor/pipeline/parsets/image_sector_pipeline.cwl
@@ -394,6 +394,14 @@ outputs:
     outputSource:
       - find_diagnostics/diagnostics
     type: File
+  - id: sector_offsets
+    outputSource:
+      - find_diagnostics/offsets
+    type: File
+  - id: sector_diagnostic_plots
+    outputSource:
+      - find_diagnostics/plots
+    type: File[]
   - id: visibilities
     outputSource:
       - prepare_imaging_data/msimg
@@ -819,5 +827,13 @@ steps:
         source: obs_filename
       - id: diagnostics_file
         source: filter/diagnostics
+      - id: facet_region_file
+{% if use_facets %}
+        source: make_region_file/region_file
+{% else %}
+        valueFrom: 'none'
+{% endif %}
     out:
       - id: diagnostics
+      - id: offsets
+      - id: plots
diff --git a/rapthor/pipeline/steps/calculate_image_diagnostics.cwl b/rapthor/pipeline/steps/calculate_image_diagnostics.cwl
index 481a78c1d3f9f010c4071c3ceb1e5ec41efa3d50..64905f1b73ceee4fc19da95d10b4baa1e8b4ef79 100644
--- a/rapthor/pipeline/steps/calculate_image_diagnostics.cwl
+++ b/rapthor/pipeline/steps/calculate_image_diagnostics.cwl
@@ -78,6 +78,18 @@ inputs:
     type: string
     inputBinding:
       position: 9
+  - id: facet_region_file
+    label: Input ds9 region file
+    doc: |
+      The filename of the input ds9 region file that defines the facets. Note
+      that when this file is unavailable, the filename can be set to a dummy
+      string, in which case it is then ignored by the script
+    type:
+      - string?
+      - File?
+    inputBinding:
+      prefix: --facet_region_file=
+      separate: false
 
 outputs:
   - id: diagnostics
@@ -87,6 +99,20 @@ outputs:
     type: File
     outputBinding:
       glob: '$(inputs.output_root).image_diagnostics.json'
+  - id: offsets
+    label: Astrometry offsets
+    doc: |
+      The astrometry offsets in RA and Dec, per facet.
+    type: File
+    outputBinding:
+      glob: '$(inputs.output_root).astrometry_offsets.json'
+  - id: plots
+    label: Diagnostic plots
+    doc: |
+      Various diagnostic plots of the photometry and astrometry.
+    type: File[]
+    outputBinding:
+      glob: '*.pdf'
 
 
 hints:
diff --git a/rapthor/scripts/calculate_image_diagnostics.py b/rapthor/scripts/calculate_image_diagnostics.py
index 758fbf1da9c62a769b89247029c777ccc1968632..826a285f9a0329f17246db883aa8e1cffcfd128e 100755
--- a/rapthor/scripts/calculate_image_diagnostics.py
+++ b/rapthor/scripts/calculate_image_diagnostics.py
@@ -7,6 +7,7 @@ from argparse import RawTextHelpFormatter
 import lsmtool
 import numpy as np
 from rapthor.lib import miscellaneous as misc
+from rapthor.lib.facet import SquareFacet, read_ds9_region_file, make_wcs, radec2xy
 from rapthor.lib.fitsimage import FITSImage
 from rapthor.lib.observation import Observation
 import casacore.tables as pt
@@ -14,9 +15,16 @@ from astropy.utils import iers
 from astropy.table import Table
 import astropy.units as u
 from astropy.coordinates import SkyCoord
-from astropy.coordinates import match_coordinates_sky
-from astropy.stats import sigma_clip
 import json
+from astropy.visualization.wcsaxes import WCSAxes
+import tempfile
+import matplotlib
+if matplotlib.get_backend() != 'Agg':
+    matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import os
+import glob
+import shutil
 
 
 # Turn off astropy's IERS downloads to fix problems in cases where compute
@@ -24,9 +32,335 @@ import json
 iers.conf.auto_download = False
 
 
+def plot_astrometry_offsets(facets, field_ra, field_dec, output_file, plot_labels=False):
+    """
+    Plots the astrometry offsets across the field
+
+    Note: The arrows indicate the direction and magnitude of the
+    corrections that should be applied to the LOFAR positions in
+    each facet to obtain a mean offset of zero
+
+    Parameters
+    ----------
+    facets : list
+        List of Facet objects
+    field_ra : float
+        RA in degrees of the field
+    field_dec : float
+        Dec in degrees of the field
+    output_file : str
+        Filename for the output plot
+    plot_labels : bool, optional
+        If True, plot the facet labels
+    """
+    wcs = make_wcs(field_ra, field_dec)
+    ra_offsets = []
+    dec_offsets = []
+    facet_patches = []
+    facet_ra = []
+    facet_dec = []
+    facet_names = []
+    for facet in facets:
+        # Note: the offsets are calculated as (LOFAR model value) - (comparison model
+        # value); e.g., a positive Dec offset indicates that the LOFAR sources
+        # lie on average to the North of the comparison source positions. Therefore, we
+        # multiply by -1 to obtain the correct arrow directions in the quiver plot
+        if facet.astrometry_diagnostics:
+            ra_offsets.append(-1 * facet.astrometry_diagnostics['meanClippedRAOffsetDeg'])
+            dec_offsets.append(-1 * facet.astrometry_diagnostics['meanClippedDecOffsetDeg'])
+            facet_ra.append(facet.ra)
+            facet_dec.append(facet.dec)
+            facet_names.append(facet.name)
+            facet_patches.append(facet.get_matplotlib_patch(wcs=wcs))
+
+    # Set up the figure. We use various values that should produce a reasonably
+    # sized figure with readable labels in most cases
+    fig = plt.figure(1, figsize=(7.66, 7))  # increase x size to ensure room for Dec label
+    plt.clf()
+    ax = WCSAxes(fig, [0.16, 0.1, 0.8, 0.8], wcs=wcs)  # set start x to give room for Dec label
+    fig.add_axes(ax)
+    RAAxis = ax.coords['ra']
+    RAAxis.set_axislabel('RA', minpad=0.75)
+    RAAxis.set_major_formatter('hh:mm:ss')
+    DecAxis = ax.coords['dec']
+    DecAxis.set_axislabel('Dec', minpad=0.75)
+    DecAxis.set_major_formatter('dd:mm:ss')
+    ax.coords.grid(color='black', alpha=0.5, linestyle='solid')
+    ax.set_title('Positional Offsets (arrows indicate direction and magnitude of correction)')
+
+    # Plot the facet polygons
+    x, y = radec2xy(wcs, facet_ra, facet_dec)
+    for i, patch in enumerate(facet_patches):
+        ax.add_patch(patch)
+        if plot_labels:
+            ax.annotate(facet_names[i], (x[i], y[i]), va='top', ha='center')
+
+    # Plot the offsets. The arrows indicate the direction and magnitude of the
+    # correction that should be applied to the LOFAR positions
+    x_offsets = np.array(ra_offsets) / wcs.wcs.cdelt[0]
+    y_offsets = np.array(dec_offsets) / wcs.wcs.cdelt[1]
+    quiver_plot = plt.quiver(x, y, x_offsets, y_offsets, units='xy', angles='xy',
+                             scale=2/3600, color='blue')
+    plt.quiverkey(quiver_plot, 0.1, 0.95, 1/3600/wcs.wcs.cdelt[0], '1 arcsec',
+                  labelpos='N', coordinates='figure')
+    plt.savefig(output_file, format='pdf')
+
+
+def fits_to_makesourcedb(catalog, reference_freq, flux_colname='Isl_Total_flux'):
+    """
+    Converts a PyBDSF catalog to a makesourcedb sky model
+
+    Note: the resulting makesourcedb catalog is a minimal catalog suitable for
+    use in photometry and astrometry comparisons and should not be used for
+    calibration
+
+    Parameters
+    ----------
+    catalog : astropy Table object
+        Input PyBDSF catalog
+    reference_freq : float
+        The reference frequency in Hz for the input catalog at which the flux
+        densities were measured
+    flux_colname : str, optional
+        The name of the column in the input catalog that contains the flux
+        density values
+
+    Returns
+    -------
+    skymodel : LSMTool skymodel object
+        The makesourcedb sky model
+    """
+    # Convert the result to makesourcedb format and write to a tempfile
+    out_lines = [f'FORMAT = Name, Type, Ra, Dec, I, ReferenceFrequency={reference_freq}\n']
+    for (name, ra, dec, flux) in zip(catalog['Source_id'], catalog['RA'],
+                                     catalog['DEC'], catalog[flux_colname]):
+        out_lines.append(f'{name}, POINT, {ra}, {dec}, {flux}, \n')
+
+    skymodel_file = tempfile.NamedTemporaryFile()
+    with open(skymodel_file.name, 'w') as f:
+        f.writelines(out_lines)
+    skymodel = lsmtool.load(skymodel_file.name)
+
+    return skymodel
+
+
+def check_photometry(obs, input_catalog, freq, min_number, comparison_skymodel=None):
+    """
+    Calculate and plot various photometry diagnostics
+
+    Parameters
+    ----------
+    obs : Observation object
+        Representative observation, used to derive pointing, etc.
+    input_catalog : str
+        Filename of the input PyBDSF FITS catalog derived from the LOFAR image
+    freq : float
+        Frequency in Hz of the LOFAR image
+    min_number : int
+        Minimum number of matched sources required for the comparisons
+    comparison_skymodel : str, optional
+        Filename of the sky model to use for the photometry (flux scale)
+        comparison (in makesourcedb format). If not given, a TGSS model is
+        downloaded
+
+    Returns
+    -------
+    photometry_diagnostics : dict
+        Photometry diagnositcs
+    """
+    # Load photometry comparison model
+    if comparison_skymodel:
+        try:
+            s_comp_photometry = lsmtool.load(comparison_skymodel)
+            name_comp_photometry = 'User supplied catalog'
+        except OSError as e:
+            # Comparison catalog not loaded successfully
+            s_comp_photometry = None
+            print('Comparison sky model could not be loaded. Error was: {}. Trying default '
+                  'sky model instead...'.format(e))
+    else:
+        s_comp_photometry = None
+    if s_comp_photometry is None:
+        try:
+            # Download a TGSS sky model around the phase center, using a 5-deg radius
+            # to ensure the field is fully covered
+            s_comp_photometry = lsmtool.load('tgss', VOPosition=[obs.ra, obs.dec], VORadius=5.0)
+            name_comp_photometry = 'TGSS'
+        except OSError as e:
+            # Comparison catalog not loaded successfully
+            print('Comparison sky model could not be loaded. Error was: {}. Skipping photometry '
+                  'check...'.format(e))
+            s_comp_photometry = None
+            photometry_diagnostics = None
+
+    # Do the photometry check
+    if s_comp_photometry:
+        # Filter the input PyBDSF FITS catalog as needed for the photometry check
+        # Sources are filtered to keep only those that:
+        #   - lie within the FWHM of the primary beam (to exclude sources with
+        #     uncertain primary beam corrections)
+        #   - have deconvolved major axis < 10 arcsec (to exclude extended sources
+        #     that may be poorly modeled)
+        catalog = Table.read(input_catalog, format='fits')
+        phase_center = SkyCoord(ra=obs.ra*u.degree, dec=obs.dec*u.degree)
+        coords_comp = SkyCoord(ra=catalog['RA'], dec=catalog['DEC'])
+        separation = phase_center.separation(coords_comp)
+        sec_el = 1.0 / np.sin(obs.mean_el_rad)
+        fwhm_deg = 1.1 * ((3.0e8 / freq) / obs.diam) * 180 / np.pi * sec_el
+        catalog = catalog[separation < fwhm_deg/2*u.degree]
+        major_axis = catalog['DC_Maj']  # degrees
+        catalog = catalog[major_axis < 10/3600]
+
+        # Convert the filtered catalog to a minimal sky model for use with LSMTool
+        # and do the comparison
+        s_pybdsf = fits_to_makesourcedb(catalog, freq)
+        s_comp_photometry.group('every')
+        if len(s_pybdsf) >= min_number:
+            photometry_diagnostics = s_pybdsf.compare(s_comp_photometry, radius='5 arcsec',
+                                                      excludeMultiple=True, make_plots=True,
+                                                      name1='LOFAR', name2=name_comp_photometry)
+        else:
+            photometry_diagnostics = None
+            print(f'Fewer than {min_number} sources found in the LOFAR image meet '
+                  'the photometry cuts (major axis < 10" and located inside the FWHM '
+                  'of the primary beam"). Skipping photometry check...')
+
+    return photometry_diagnostics
+
+
+def check_astrometry(obs, input_catalog, image, facet_region_file, min_number,
+                     output_root, comparison_skymodel=None):
+    """
+    Calculate and plot various astrometry diagnostics
+
+    Parameters
+    ----------
+    obs : Observation object
+        Representative observation, used to derive pointing, etc.
+    input_catalog : str
+        Filename of the input PyBDSF FITS catalog derived from the LOFAR image
+    image : FITSImage object
+        The LOFAR image, used to derive frequency, coverage, etc.
+    facet_region_file : str, optional
+        Filename of the facet region file (in ds9 format) that defines the
+        facets used in imaging
+    min_number : int
+        Minimum number of matched sources required for the comparisons
+    output_root : str
+        Root of the filename for the output files
+    comparison_skymodel : str, optional
+        Filename of the sky model to use for the photometry (flux scale)
+        comparison (in makesourcedb format). If not given, a Pan-STARRS model is
+        downloaded
+
+    Returns
+    -------
+    astrometry_diagnostics : dict
+        Astrometry diagnositcs
+    """
+    # Load and filter the input PyBDSF FITS catalog as needed for the astrometry check
+    # Sources are filtered to keep only those that:
+    #   - have deconvolved major axis < 10 arcsec (to exclude extended sources
+    #     that may be poorly modeled)
+    #   - have errors on RA and Dec of < 2 arcsec (to exclude sources
+    #     with high positional uncertainties)
+    catalog = Table.read(input_catalog, format='fits')
+    major_axis = catalog['DC_Maj']  # degrees
+    catalog = catalog[major_axis < 10/3600]
+    ra_error = catalog['E_RA']  # degrees
+    catalog = catalog[ra_error < 2/3600]
+    dec_error = catalog['E_DEC']  # degrees
+    catalog = catalog[dec_error < 2/3600]
+
+    # Do the astrometry check
+    if len(catalog) >= min_number:
+        max_search_cone_radius = 0.5  # deg; Pan-STARRS search limit
+        if facet_region_file is not None and os.path.isfile(facet_region_file):
+            facets = read_ds9_region_file(facet_region_file)
+        else:
+            # Use a single rectangular facet centered on the phase center
+            ra = obs.ra
+            dec = obs.dec
+            image_width = max(image.img_data.shape[-2:]) * abs(image.img_hdr['CDELT1'])
+            width = min(max_search_cone_radius*2, image_width)
+            facets = [SquareFacet('field', ra, dec, width)]
+
+        # Convert the filtered catalog into a minimal sky model for use with LSMTool
+        s_pybdsf = fits_to_makesourcedb(catalog, image.freq)
+
+        # Loop over the facets, performing the astrometry checks for each
+        if comparison_skymodel:
+            try:
+                s_comp_astrometry = lsmtool.load(comparison_skymodel)
+                s_comp_astrometry.group('every')
+            except OSError as e:
+                # Comparison catalog not loaded successfully
+                s_comp_astrometry = None
+                print('Comparison sky model could not be loaded. Error was: {}. Trying default '
+                      'sky model instead...'.format(e))
+        else:
+            s_comp_astrometry = None
+
+        astrometry_diagnostics = {'facet_name': [],
+                                  'meanRAOffsetDeg': [],
+                                  'stdRAOffsetDeg': [],
+                                  'meanClippedRAOffsetDeg': [],
+                                  'stdClippedRAOffsetDeg': [],
+                                  'meanDecOffsetDeg': [],
+                                  'stdDecOffsetDeg': [],
+                                  'meanClippedDecOffsetDeg': [],
+                                  'stdClippedDecOffsetDeg': []}
+        for facet in facets:
+            facet.set_skymodel(s_pybdsf.copy())
+            facet.find_astrometry_offsets(s_comp_astrometry, min_number=min_number)
+            if len(facet.astrometry_diagnostics):
+                astrometry_diagnostics['facet_name'].append(facet.name)
+                astrometry_diagnostics['meanRAOffsetDeg'].append(facet.astrometry_diagnostics['meanRAOffsetDeg'])
+                astrometry_diagnostics['stdRAOffsetDeg'].append(facet.astrometry_diagnostics['stdRAOffsetDeg'])
+                astrometry_diagnostics['meanClippedRAOffsetDeg'].append(facet.astrometry_diagnostics['meanClippedRAOffsetDeg'])
+                astrometry_diagnostics['stdClippedRAOffsetDeg'].append(facet.astrometry_diagnostics['stdClippedRAOffsetDeg'])
+                astrometry_diagnostics['meanDecOffsetDeg'].append(facet.astrometry_diagnostics['meanDecOffsetDeg'])
+                astrometry_diagnostics['stdDecOffsetDeg'].append(facet.astrometry_diagnostics['stdDecOffsetDeg'])
+                astrometry_diagnostics['meanClippedDecOffsetDeg'].append(facet.astrometry_diagnostics['meanClippedDecOffsetDeg'])
+                astrometry_diagnostics['stdClippedDecOffsetDeg'].append(facet.astrometry_diagnostics['stdClippedDecOffsetDeg'])
+
+        # Save and plot the per-facet offsets
+        if len(astrometry_diagnostics['facet_name']):
+            with open(output_root+'.astrometry_offsets.json', 'w') as fp:
+                json.dump(astrometry_diagnostics, fp)
+            ra = obs.ra
+            dec = obs.dec
+            plot_astrometry_offsets(facets, ra, dec, output_root+'.astrometry_offsets.pdf')
+
+            # Calculate mean offsets
+            mean_astrometry_diagnostics = {'meanRAOffsetDeg': np.mean(astrometry_diagnostics['meanRAOffsetDeg']),
+                                           'stdRAOffsetDeg': np.mean(astrometry_diagnostics['stdRAOffsetDeg']),
+                                           'meanClippedRAOffsetDeg': np.mean(astrometry_diagnostics['meanClippedRAOffsetDeg']),
+                                           'stdClippedRAOffsetDeg': np.mean(astrometry_diagnostics['stdClippedRAOffsetDeg']),
+                                           'meanDecOffsetDeg': np.mean(astrometry_diagnostics['meanDecOffsetDeg']),
+                                           'stdDecOffsetDeg': np.mean(astrometry_diagnostics['stdDecOffsetDeg']),
+                                           'meanClippedDecOffsetDeg': np.mean(astrometry_diagnostics['meanClippedDecOffsetDeg']),
+                                           'stdClippedDecOffsetDeg': np.mean(astrometry_diagnostics['stdClippedDecOffsetDeg'])}
+        else:
+            # Write dummy files
+            mean_astrometry_diagnostics = None
+            with open(output_root+'.astrometry_offsets.json', 'w') as fp:
+                fp.writelines('Astrometry diagnostics could not be determined. Please see the logs for details')
+            with open(output_root+'.astrometry_offsets.pdf', 'w') as fp:
+                fp.writelines('Astrometry diagnostics could not be determined. Please see the logs for details')
+    else:
+        mean_astrometry_diagnostics = None
+        print(f'Fewer than {min_number} sources found in the LOFAR image meet the '
+              'astrometry cuts (major axis < 10" with positional errors < 2"). '
+              'Skipping the astromety check...')
+
+    return mean_astrometry_diagnostics
+
+
 def main(flat_noise_image, flat_noise_rms_image, true_sky_image, true_sky_rms_image,
          input_catalog, input_skymodel, obs_ms, diagnostics_file, output_root,
-         comparison_skymodel=None):
+         facet_region_file=None, photometry_comparison_skymodel=None,
+         astrometry_comparison_skymodel=None, min_number=5):
     """
     Calculate various image diagnostics
 
@@ -52,9 +386,18 @@ def main(flat_noise_image, flat_noise_rms_image, true_sky_image, true_sky_rms_im
         by the sky model filtering script
     output_root : str
         Root of the filename for the output files
-    comparison_skymodel : str, optional
-        Filename of the sky model to use for flux scale and astrometry
-        comparisons. If not given, a TGSS model is downloaded
+    facet_region_file : str, optional
+        Filename of the facet region file (in ds9 format) that defines the
+        facets used in imaging
+    photometry_comparison_skymodel : str, optional
+        Filename of the sky model to use for the photometry (flux scale)
+        comparison (in makesourcedb format). If not given, a TGSS model is
+        downloaded
+    astrometry_comparison_skymodel : str, optional
+        Filename of the sky model to use for the astrometry comparison (in
+        makesourcedb format). If not given, a Pan-STARRS model is downloaded
+    min_number : int, optional
+        Minimum number of matched sources required for the comparisons
     """
     # Select the best MS
     if isinstance(obs_ms, str):
@@ -104,114 +447,35 @@ def main(flat_noise_image, flat_noise_rms_image, true_sky_image, true_sky_rms_im
                        'freq': img_true_sky.freq,
                        'beam_fwhm': img_true_sky.beam})
 
-    # Load input sky model and comparison model
-    s_in = lsmtool.load(input_skymodel)
-    if comparison_skymodel is None:
-        # Download a TGSS sky model around the midpoint of the input sky model,
-        # using a 5-deg radius to ensure the field is fully covered
-        _, _, midRA, midDec = s_in._getXY()
-        try:
-            s_comp = lsmtool.load('tgss', VOPosition=[midRA, midDec], VORadius=5.0)
-        except OSError:
-            # Comparison catalog not downloaded successfully
-            s_comp = None
-    else:
-        s_comp = lsmtool.load(comparison_skymodel)
-
-    # Do cross matching of the comparison model with the PyBDSF-derived one
-    if s_comp and s_in:
-        # Write the comparison catalog to FITS file for use with astropy
-        catalog_comp_filename = output_root+'.comparison.fits'
-        s_comp.table.columns['Ra'].format = None
-        s_comp.table.columns['Dec'].format = None
-        s_comp.table.columns['I'].format = None
-        s_comp.write(catalog_comp_filename, format='fits', clobber=True)
-        if 'ReferenceFrequency' in s_comp.getColNames():
-            freq_comp = np.mean(s_comp.getColValues('ReferenceFrequency'))
-        else:
-            freq_comp = s_comp.table.meta['ReferenceFrequency']
+    # Do the photometry check and update the ouput dict
+    result = check_photometry(obs_list[beam_ind], input_catalog, img_true_sky.freq,
+                              min_number, comparison_skymodel=photometry_comparison_skymodel)
+    if result is not None:
+        cwl_output.update(result)
 
-        # Read in PyBDSF-derived and comparison catalogs and filter to keep
-        # only comparison sources within a radius of FWHM / 2 of phase center
-        catalog = Table.read(input_catalog, format='fits')
-        catalog_comp = Table.read(catalog_comp_filename, format='fits')
-        obs = obs_list[beam_ind]
-        phase_center = SkyCoord(ra=obs.ra*u.degree, dec=obs.dec*u.degree)
-        coords_comp = SkyCoord(ra=catalog_comp['Ra'], dec=catalog_comp['Dec'])
-        separation = phase_center.separation(coords_comp)
-        sec_el = 1.0 / np.sin(obs.mean_el_rad)
-        fwhm_deg = 1.1 * ((3.0e8 / img_true_sky.freq) / obs.diam) * 180 / np.pi * sec_el
-        catalog_comp = catalog_comp[separation < fwhm_deg/2*u.degree]
-
-        # Cross match the coordinates and keep only matches that have a
-        # separation of 5 arcsec or less
-        #
-        # Note: we require at least 10 sources for the comparison, as using
-        # fewer can give unreliable estimates
-        min_number = 10
-        if len(catalog_comp) >= min_number:
-            coords = SkyCoord(ra=catalog['RA'], dec=catalog['DEC'])
-            coords_comp = SkyCoord(ra=catalog_comp['Ra'], dec=catalog_comp['Dec'])
-            idx, sep2d, _ = match_coordinates_sky(coords_comp, coords)
-            constraint = sep2d < 5*u.arcsec
-            catalog_comp = catalog_comp[constraint]
-            catalog = catalog[idx[constraint]]
-
-            # Find the mean flux ratio (input / comparison). We use the total island
-            # flux density from PyBDSF as it gives less scatter than the Gaussian
-            # fluxes when comparing to a lower-resolution catalog (such as TGSS).
-            # Lastly, a correction factor is used to correct for the potentially
-            # different frequencies of the LOFAR image and the comparison catalog
-            if len(catalog_comp) >= min_number:
-                alpha = -0.7  # adopt a typical spectral index
-                freq_factor = (freq_comp / img_true_sky.freq)**alpha
-                ratio = catalog['Isl_Total_flux'] / catalog_comp['I'] * freq_factor
-                meanRatio = np.mean(ratio)
-                stdRatio = np.std(ratio)
-                clippedRatio = sigma_clip(ratio)
-                meanClippedRatio = np.mean(clippedRatio)
-                stdClippedRatio = np.std(clippedRatio)
-                stats = {'meanRatio_pybdsf': meanRatio,
-                         'stdRatio_pybdsf': stdRatio,
-                         'meanClippedRatio_pybdsf': meanClippedRatio,
-                         'stdClippedRatio_pybdsf': stdClippedRatio}
-                cwl_output.update(stats)
-
-        # Group the comparison sky model into sources using a FWHM of 40 arcsec,
-        # the approximate TGSS resolution
-        s_comp.group('threshold', FWHM='40.0 arcsec', threshold=0.05)
-
-        # Keep POINT-only sources
-        source_type = s_comp.getColValues('Type')
-        patch_names = s_comp.getColValues('Patch')
-        non_point_patch_names = set(patch_names[np.where(source_type != 'POINT')])
-        ind = []
-        for patch_name in non_point_patch_names:
-            ind.extend(s_comp.getRowIndex(patch_name))
-        s_comp.remove(np.array(ind))
-
-        # Check if there is a sufficient number of sources to do the comparison with.
-        # If there is, do it and append the resulting diagnostics dict to the
-        # existing one
-        #
-        # Note: the various ratios are all calculated as (s_in / s_comp) and the
-        # differences as (s_in - s_comp). If there are no successful matches,
-        # the compare() method returns None
-        if (s_comp
-                and len(s_in.getPatchNames()) >= min_number
-                and len(s_comp.getPatchNames()) >= min_number):
-            flux_astrometry_diagnostics = s_in.compare(s_comp, radius='5 arcsec',
-                                                       excludeMultiple=True, make_plots=False)
-            if flux_astrometry_diagnostics is not None:
-                cwl_output.update(flux_astrometry_diagnostics)
-
-    # Write out the diagnostics
+    # Do the astrometry check and update the ouput dict
+    result = check_astrometry(obs_list[beam_ind], input_catalog, img_true_sky, facet_region_file,
+                              min_number, output_root, comparison_skymodel=astrometry_comparison_skymodel)
+    if result is not None:
+        cwl_output.update(result)
+
+    # Write out the full diagnostics
     with open(output_root+'.image_diagnostics.json', 'w') as fp:
         json.dump(cwl_output, fp)
 
+    # Adjust plot filenames as needed to include output_root (those generated by LSMTool
+    # lack it)
+    diagnostic_plots = glob.glob(os.path.join('.', '*.pdf'))
+    for src_filename in diagnostic_plots:
+        if not os.path.basename(src_filename).startswith(output_root):
+            dst_filename = os.path.join('.', f'{output_root}.' + os.path.basename(src_filename))
+            if os.path.exists(dst_filename):
+                os.remove(dst_filename)
+            shutil.copy(src_filename, dst_filename)
+
 
 if __name__ == '__main__':
-    descriptiontext = "Analyze source catalogs.\n"
+    descriptiontext = "Calculate image photometry and astrometry diagnostics.\n"
 
     parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
     parser.add_argument('flat_noise_image', help='Filename of flat-noise FITS image')
@@ -223,7 +487,13 @@ if __name__ == '__main__':
     parser.add_argument('obs_ms', help='Filename of observation MS')
     parser.add_argument('diagnostics_file', help='Filename of diagnostics JSON file')
     parser.add_argument('output_root', help='Root of output files')
+    parser.add_argument('--facet_region_file', help='Filename of ds9 facet region file', type=str, default=None)
+    parser.add_argument('--photometry_comparison_skymodel', help='Filename of photometry sky model', type=str, default=None)
+    parser.add_argument('--astrometry_comparison_skymodel', help='Filename of astrometry sky model', type=str, default=None)
+    parser.add_argument('--min_number', help='Minimum number of sources for diagnostics', type=int, default=5)
 
     args = parser.parse_args()
     main(args.flat_noise_image, args.flat_noise_rms_image, args.true_sky_image, args.true_sky_rms_image,
-         args.input_catalog, args.input_skymodel, args.obs_ms, args.diagnostics_file, args.output_root)
+         args.input_catalog, args.input_skymodel, args.obs_ms, args.diagnostics_file, args.output_root,
+         facet_region_file=args.facet_region_file, photometry_comparison_skymodel=args.photometry_comparison_skymodel,
+         astrometry_comparison_skymodel=args.astrometry_comparison_skymodel, min_number=args.min_number)
diff --git a/requirements.txt b/requirements.txt
index 4d213b016ca747316f3f72887330c8ab53e17eff..82548987ded65c9de8acdaff05f7f296705b6f12 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,6 +10,7 @@ psutil>=5.6.6
 python-casacore
 python-dateutil
 reproject
+requests
 Rtree
 scipy
 shapely