Ameliorations WCS

Alain Klotz
1 parent ad2b1e4a
Showing 5 changed files with 180 additions and 9 deletions Show diff stats
.gitignore
resources/solar_system/de421.bsp
src/guitastro/.gitignore
src/guitastro/ephemeris.py
src/guitastro/wcs.py
 __pycache__
-de421.bsp
 photometrie1.txt
 requirements.txt
 build
@@ -0,0 +1 @@
+de421.bsp
@@ -1980,10 +1980,10 @@ if __name__ == &quot;__main__&quot;:
  
     if example == 16:
         """
-        Load a matrix file generated by PyROS
+        Load a matrix file schedule generated by PyROS
         """
         eph = Ephemeris()
-        fname = r"C:\srv\develop\pyros_soft\pyros\data\sequence\P001\20230831\scheduler_schedule.txt"
+        fname = r"C:\srv\develop\pyros_soft\DATA\PRODUCTS\TESTS\LONGSTORE\sequences\P001\20230925\scheduler_schedule.txt"
         eph.skdl_plot_schedule(fname)
  
  
@@ -800,6 +800,134 @@ class Wcs(WcsException, GuitastroTools):
         self.wcs2skyfield(naxis1, naxis2)
         return self._skyfield;
  
+    def skyfield2fov(self, skyfield:dict) ->tuple:
+        """ Compute the image field of view profected on R.A., Dec. from a skyfield.
+
+        Args:
+            skyfield: Input values.
+
+        Returns:
+            Two elements defining the field of view:
+
+                * fovra: (float) Field of view along Right Ascension (unit degrees)
+                * fovdec: (float) Field of view along Declination (unit degrees)
+        """
+        cd11 = math.degrees(skyfield["cd11rad"])
+        cd12 = math.degrees(skyfield["cd12rad"])
+        cd21 = math.degrees(skyfield["cd21rad"])
+        cd22 = math.degrees(skyfield["cd22rad"])
+        naxis1 = skyfield["naxis1"]
+        naxis2 = skyfield["naxis2"]
+        cdelt1, cdelt2, crota2 = self.cd2cdelt(cd11, cd12, cd21, cd22)
+        crota2r = math.radians(crota2)
+        cosrota = math.cos(crota2r)
+        sinrota = math.sin(crota2r)
+        fov1 = cdelt1*naxis1
+        fov2 = cdelt2*naxis2
+        fovra  = abs(fov1*cosrota) + abs(fov2*sinrota)
+        fovdec = abs(fov1*sinrota) + abs(fov2*cosrota)
+        return fovra, fovdec
+
+    def is_skyfield_in_radecfov(self, ra:float, dec:float, fovra:float, fovdec:float, coscrval1:float, sincrval1: float, coscrval2:float, sincrval2: float) ->bool:
+        """Check if (ra,dec) target is perhaps inside a skyfields.
+
+        Skyfield is defined by 12 elements.
+        Four of them are required as input of this method: coscrval1, sincrval1, coscrval2, sincrval2.
+        (ra, dec) target must be related to the same equinox than crval1, crval2 coordinates.
+        The code of this method can be adapted in programs that manage a SQL database.
+        The algorithm is rapid but not rigourous.
+        The goal is to exclude skyfields well outside the field of view.
+        Before using is_skyfield_in_radecfov() use skyfield2fov().
+        After using is_skyfield_in_radecfov() use skyfields2indexes().
+
+        Args:
+            ra: Right ascension (degrees) of the target
+            dec: Declination (degrees) of the target
+            fov: Field of view (degrees) of the images stored in skyfield
+            coscrval1: cosinus of the CRVAL1 coordinate (degrees) of the skyfield
+            sincrval1: sinus of the CRVAL1 coordinate (degrees) of the skyfield
+            coscrval2: cosinus of the CRVAL2 coordinate (degrees) of the skyfield
+            sincrval2: sinus of the CRVAL2 coordinate (degrees) of the skyfield
+
+        Returns:
+            Bolean value is True if the target should be inside the field of view.
+        """
+        decr = math.radians(dec)
+        fovracos = fovra/math.cos(decr)
+        if fovracos > 180:
+            fovracos = 180
+        # --- ra
+        rar = math.radians(ra+fovracos/2)
+        cosrap= math.cos(rar)
+        sinrap= math.sin(rar)
+        rar = math.radians(ra-fovracos/2)
+        cosram= math.cos(rar)
+        sinram= math.sin(rar)
+        if cosrap < cosram:
+            cosrap , cosram = cosram , cosrap
+        if sinrap < sinram:
+            sinrap , sinram = sinram , sinrap
+        if sinram < 0 and sinrap >= 0:
+            if cosram > 0:
+                cosrap = 1
+            else:
+                cosrap = -1
+        if cosram < 0 and cosrap >= 0:
+            if sinram > 0:
+                sinrap = 1
+            else:
+                sinrap = -1
+        # --- dec
+        decr = math.radians(dec+fovdec/2)
+        cosdecp= math.cos(decr)
+        sindecp= math.sin(decr)
+        decr = math.radians(dec-fovdec/2)
+        cosdecm= math.cos(decr)
+        sindecm= math.sin(decr)
+        # --- re-arange min, max
+        if cosdecp < cosdecm:
+            cosdecp , cosdecm = cosdecm , cosdecp
+        if sindecp < sindecm:
+            sindecp , sindecm = sindecm , sindecp
+        if sindecm < 0 and sindecp >= 0:
+            if cosdecm > 0:
+                cosdecp = 1
+            else:
+                cosdecp = -1
+        if cosdecm < 0 and cosdecp >= 0:
+            if sindecm > 0:
+                sindecp = 1
+            else:
+                sindecp = -1
+        # --- SQL request
+        # sql_query = "WHERE"
+        # sql_query += f" `coscrval1` >= {cosram}"
+        # sql_query += f" AND `coscrval1` <= {cosrap}"
+        # sql_query += f" AND `sincrval1` >= {sinram}"
+        # sql_query += f" AND `sincrval1` <= {sinrap}"
+        # sql_query += f" AND `coscrval2` >= {cosdecm}"
+        # sql_query += f" AND `coscrval2` <= {cosdecp}"
+        # sql_query += f" AND `sincrval2` >= {sindecm}"
+        # sql_query += f" AND `sincrval2` <= {sindecp}"
+        # --- test is true only if 4 crval* are inside the limits
+        if coscrval1 < cosram:
+            return False
+        if coscrval1 > cosrap:
+            return False
+        if sincrval1 < sinram:
+            return False
+        if sincrval1 > sinrap:
+            return False
+        if coscrval2 < cosdecm:
+            return False
+        if coscrval2 > cosdecp:
+            return False
+        if sincrval2 < sindecm:
+            return False
+        if sincrval2 > sindecp:
+            return False
+        return True
+
     def skyfields2indexes(self, ra:float, dec:float, skyfields:list):
         """Check if (ra,dec) is inside a list of skyfields.
  
@@ -846,7 +974,7 @@ class Wcs(WcsException, GuitastroTools):
             cd22 = skyfield["cd22rad"]
             crp1 = skyfield["crpix1"]
             crp2 = skyfield["crpix2"]
-            naxis1 = skyfield["naxis2"]
+            naxis1 = skyfield["naxis1"]
             naxis2 = skyfield["naxis2"]
             ccca = cc*ca # cos(dec)*cos(ra) * cos(ra0)
             cssa = cs*sa # cos(dec)*sin(ra) * sin(ra0)
@@ -878,8 +1006,8 @@ class Wcs(WcsException, GuitastroTools):
  
 if __name__ == "__main__":
  
-    default = 3
-    example = input(f"Select the example (0 to 3) ({default}) ")
+    default = 4
+    example = input(f"Select the example (0 to 4) ({default}) ")
     try:
         example = int(example)
     except:
@@ -912,14 +1040,57 @@ if __name__ == &quot;__main__&quot;:
  
     if example == 3:
         """
-        Use of skyfield
+        Use of skyfields2indexes
         """
         skyfields = []
         wcs = Wcs()
-        wcs.optic(512, 512, "6h00m", "+30d", 9, 9, 1, 0)
+        wcs.optic(512, 512, "8h00m", "+34d", 9, 9, 1, 0)
         skyfield = wcs.skyfield_get(1024, 1024)
         skyfields.append(skyfield)
-        wcs.optic(512, 512, "8h00m", "+34d", 9, 9, 1, 0)
+        wcs.optic(512, 512, "6h00m", "+30d", 9, 9, 1, 0)
         skyfield = wcs.skyfield_get(1024, 1024)
         skyfields.append(skyfield)
         indexes = wcs.skyfields2indexes(90.2, 30.04, skyfields)
+        for index in indexes:
+            i, x, y = index
+            ra, dec = wcs.xy2radec(x, y)
+            print(f"Target {i} position in the image {x=:.2f} {y=:.2f}. R.A.={ra:.2f} Dec.={dec:.2f}")
+
+    if example == 4:
+        """
+        Use of is_skyfield_in_radecfov
+        """
+        wcs = Wcs()
+        # - Skyfield definition for the test
+        wcs.optic(512, 512, "0h00m", "+31d", 9, 9, 1, 0)
+        skyfield = wcs.skyfield_get(1024, 1024)
+        fovra, fovdec = wcs.skyfield2fov(skyfield)
+        print("=== Skyfield characteristics:")
+        print(f"Image field of view {fovra:.2f} x {fovdec:.2f} degrees")
+        crval1 = wcs.getval('CRVAL1')[0]
+        crval2 = wcs.getval('CRVAL2')[0]
+        print(f"Image center R.A.={crval1:.2f} Dec.={crval2:.2f} degrees")
+        # - Get the four elements
+        coscrval1 = skyfield['coscrval1']
+        sincrval1 = skyfield['sincrval1']
+        coscrval2 = skyfield['coscrval2']
+        sincrval2 = skyfield['sincrval2']
+        # - Target to check if inside the image defined by skyfield
+        print("=== Target characteristics:")
+        ra = 0.3
+        dec = 31.22
+        res = wcs.is_skyfield_in_radecfov(ra, dec, fovra, fovdec, coscrval1, sincrval1, coscrval2, sincrval2)
+        resp = f"Target R.A.={ra:.2f} Dec.={dec:.2f} is "
+        if res:
+            resp += "inside the image."
+        else:
+            resp += "outside the image."
+        print(f"{resp}")
+        if res:
+            skyfields = []
+            skyfields.append(skyfield)
+            indexes = wcs.skyfields2indexes(ra, dec, skyfields)
+            for index in indexes:
+                i, x, y = index
+                print(f"Target position in the image {x=:.2f} {y=:.2f}")
+