This is gnuastro.info, produced by makeinfo version 7.2 from
gnuastro.texi.

This book documents version 0.23.119-586f of the GNU Astronomy Utilities
(Gnuastro).  Gnuastro provides various programs and libraries for
astronomical data manipulation and analysis.

   Copyright © 2015-2025 Free Software Foundation, Inc.

     Permission is granted to copy, distribute and/or modify this
     document under the terms of the GNU Free Documentation License,
     Version 1.3 or any later version published by the Free Software
     Foundation; with no Invariant Sections, no Front-Cover Texts, and
     no Back-Cover Texts.  A copy of the license is included in the
     section entitled "GNU Free Documentation License".
INFO-DIR-SECTION Astronomy
START-INFO-DIR-ENTRY
* Gnuastro: (gnuastro).       GNU Astronomy Utilities.
* libgnuastro: (gnuastro)Gnuastro library. Full Gnuastro library doc.

* help-gnuastro: (gnuastro)help-gnuastro mailing list. Getting help.

* bug-gnuastro: (gnuastro)Report a bug. How to report bugs

* Arithmetic: (gnuastro)Arithmetic. Arithmetic operations on pixels.
* astarithmetic: (gnuastro)Invoking astarithmetic. Options to Arithmetic.

* BuildProgram: (gnuastro)BuildProgram. Compile and run programs using Gnuastro's library.
* astbuildprog: (gnuastro)Invoking astbuildprog. Options to BuildProgram.

* ConvertType: (gnuastro)ConvertType. Convert different file types.
* astconvertt: (gnuastro)Invoking astconvertt. Options to ConvertType.

* Convolve: (gnuastro)Convolve. Convolve an input file with kernel.
* astconvolve: (gnuastro)Invoking astconvolve. Options to Convolve.

* CosmicCalculator: (gnuastro)CosmicCalculator. For cosmological params.
* astcosmiccal: (gnuastro)Invoking astcosmiccal. Options to CosmicCalculator.

* Crop: (gnuastro)Crop. Crop region(s) from image(s).
* astcrop: (gnuastro)Invoking astcrop. Options to Crop.

* Fits: (gnuastro)Fits. View and manipulate FITS extensions and keywords.
* astfits: (gnuastro)Invoking astfits. Options to Fits.

* MakeCatalog: (gnuastro)MakeCatalog. Make a catalog from labeled image.
* astmkcatalog: (gnuastro)Invoking astmkcatalog. Options to MakeCatalog.

* MakeProfiles: (gnuastro)MakeProfiles. Make mock profiles.
* astmkprof: (gnuastro)Invoking astmkprof. Options to MakeProfiles.

* Match: (gnuastro)Match. Match two separate catalogs.
* astmatch: (gnuastro)Invoking astmatch. Options to Match.

* NoiseChisel: (gnuastro)NoiseChisel. Detect signal in noise.
* astnoisechisel: (gnuastro)Invoking astnoisechisel. Options to NoiseChisel.

* Segment: (gnuastro)Segment. Segment detections based on signal structure.
* astsegment: (gnuastro)Invoking astsegment. Options to Segment.

* Query: (gnuastro)Query. Access remote databases for downloading data.
* astquery: (gnuastro)Invoking astquery. Options to Query.

* Statistics: (gnuastro)Statistics. Get image Statistics.
* aststatistics: (gnuastro)Invoking aststatistics. Options to Statistics.

* Table: (gnuastro)Table. Read and write FITS binary or ASCII tables.
* asttable: (gnuastro)Invoking asttable. Options to Table.

* Warp: (gnuastro)Warp. Warp a dataset to a new grid.
* astwarp: (gnuastro)Invoking astwarp. Options to Warp.

* astscript: (gnuastro)Installed scripts. Gnuastro's installed scripts.
* astscript-ds9-region: (gnuastro)Invoking astscript-ds9-region. Options to this script
* astscript-fits-view: (gnuastro)Invoking astscript-fits-view. Options to this script
* astscript-pointing-simulate: (gnuastro)Invoking astscript-pointing-simulate. Options to this script
* astscript-psf-scale-factor: (gnuastro)Invoking astscript-psf-scale-factor. Options to this script
* astscript-psf-select-stars: (gnuastro)Invoking astscript-psf-select-stars. Options to this script
* astscript-psf-stamp: (gnuastro)Invoking astscript-psf-stamp. Options to this script
* astscript-psf-subtract: (gnuastro)Invoking astscript-psf-subtract. Options to this script
* astscript-psf-unite: (gnuastro)Invoking astscript-psf-unite. Options to this script
* astscript-radial-profile: (gnuastro)Invoking astscript-radial-profile. Options to this script
* astscript-sort-by-night: (gnuastro)Invoking astscript-sort-by-night. Options to this script
* astscript-zeropoint: (gnuastro)Invoking astscript-zeropoint. Options to this script
END-INFO-DIR-ENTRY


File: gnuastro.info,  Node: Query,  Prev: Table,  Up: Data containers

5.4 Query
=========

There are many astronomical databases available for downloading
astronomical data.  Most follow the International Virtual Observatory
Alliance (IVOA, <https://ivoa.net>) standards (and in particular the
Table Access Protocol, or TAP(1)). With TAP, it is possible to submit
your queries via a command-line downloader (for example, ‘curl’) to only
get specific tables, targets (rows in a table) or measurements (columns
in a table): you do not have to download the full table (which can be
very large in some cases)!  These customizations are done through the
Astronomical Data Query Language (ADQL(2)).

   Therefore, if you are sufficiently familiar with TAP and ADQL, you
can easily custom-download any part of an online dataset.  However, you
also need to keep a record of the URLs of each database and in many
cases, the commands will become long and hard/buggy to type on the
command-line.  On the other hand, most astronomers do not know TAP or
ADQL at all, and are forced to go to the database's web page which is
slow (it needs to download so many images, and has too much annoying
information), requires manual interaction (further making it slow and
buggy), and cannot be automated.

   Gnuastro's Query program is designed to be the middle-man in this
process: it provides a simple high-level interface to let you specify
your constraints on what you want to download.  It then internally
constructs the command to download the data based on your inputs and
runs it to download your desired data.  Query also prints the full
command before it executes it (if not called with ‘--quiet’).  Also, if
you ask for a FITS output table, the full command is written into its
0-th extension along with other input parameters to query (all Gnuastro
programs generally keep their input configuration parameters as FITS
keywords in the zero-th output).  You can see it with Gnuastro's Fits
program, like below:

     $ astfits query-output.fits -h0

   With the full command used to download the dataset, you only need a
minimal knowledge of ADQL to do lower-level customizations on your
downloaded dataset.  You can simply copy that command and change the
parts of the query string you want: ADQL is very powerful!  For example,
you can ask the server to do mathematical operations on the columns and
apply selections after those operations, or combine/match multiple
datasets.  We will try to add high-level interfaces for such
capabilities, but generally, do not limit yourself to the high-level
operations (that cannot cover everything!).

* Menu:

* Available databases::         List of available databases to Query.
* Invoking astquery::           Inputs, outputs and configuration of Query.

   ---------- Footnotes ----------

   (1) <https://ivoa.net/documents/TAP>

   (2) <https://ivoa.net/documents/ADQL>


File: gnuastro.info,  Node: Available databases,  Next: Invoking astquery,  Prev: Query,  Up: Query

5.4.1 Available databases
-------------------------

The current list of databases supported by Query are listed at the end
of this section.  To get the list of available datasets within each
database, you can use the ‘--information’ option.  for example, with the
command below you can get a list of the roughly 100 datasets that are
available within the ESA Gaia server with their description:

     $ astquery gaia --information

However, other databases like VizieR host many more datasets (tens of
thousands!).  Therefore it is very inconvenient to get the _full_
information every time you want to find your dataset of interest (the
full metadata file VizieR is more than 20Mb).  In such cases, you can
limit the downloaded and displayed information with the ‘--limitinfo’
option.  For example, with the first command below, you can get all
datasets relating to the MUSE (an instrument on the Very Large
Telescope), and those that include Roland Bacon (Principle Investigator
of MUSE) as an author (‘Bacon, R.’).  Recall that ‘-i’ is the short
format of ‘--information’.

     $ astquery vizier -i --limitinfo=MUSE
     $ astquery vizier -i --limitinfo="Bacon R."

   Once you find the recognized name of your desired dataset, you can
see the column information of that dataset with adding the dataset name.
For example, with the command below you can see the column metadata in
the ‘J/A+A/608/A2/udf10’ dataset (one of the datasets in the search
above) using this command:

     $ astquery vizier --dataset=J/A+A/608/A2/udf10 -i

   For very popular datasets of a database, Query provides an
easier-to-remember short name that you can feed to ‘--dataset’.  This
short name will map to the officially recognized name of the dataset on
the server.  In this mode, Query will also set positional columns
accordingly.  For example, most VizieR datasets have an ‘RAJ2000’ column
(the RA and the epoch of 2000) so it is the default RA column name for
coordinate search (using ‘--center’ or ‘--overlapwith’).  However, some
datasets do not have this column (for example, SDSS DR12).  So when you
use the short name and Query knows about this dataset, it will
internally set the coordinate columns that SDSS DR12 has: ‘RA_ICRS’ and
‘DEC_ICRS’.  Recall that you can always change the coordinate columns
with ‘--ccol’.

   For example, in the VizieR and Gaia databases, the recognized name
for data release 3 data is respectively ‘I/355/gaiadr3’ and
‘gaiadr3.gaia_source’.  These technical names are hard to remember.
Therefore Query provides ‘gaiadr3’ (for VizieR) and ‘dr3’ (for ESA's
Gaia database) shortcuts which you can give to ‘--dataset’ instead.
They will be internally mapped to the fully recognized name by Query.
In the list below that describes the available databases, the available
short names, that are recognized for each, are also listed.

*Not all datasets support TAP:* Large databases like VizieR have TAP
access for all their datasets.  However, smaller databases have not
implemented TAP for all their tables.  Therefore some datasets that are
searchable in their web interface may not be available for a TAP search.
To see the full list of TAP-ed datasets in a database, use the
‘--information’ (or ‘-i’) option with the dataset name like the command
below.

     $ astquery astron -i

If your desired dataset is not in this list, but has web-access, contact
the database maintainers and ask them to add TAP access for it.  After
they do it, you should see the name added to the output list of the
command above.

   The list of databases recognized by Query (and their names in Query)
is described below.  Since Query is a new member of the Gnuastro family
(first available in Gnuastro 0.14), this list will hopefully grow
significantly in the next releases.  If you have any particular datasets
in mind, please let us know by sending an email to
‘bug-gnuastro@gnu.org’.  If the dataset supports IVOA's TAP (Table
Access Protocol), it should be very easy to add.

‘astron’
     The ASTRON Virtual Observatory service (<https://vo.astron.nl>) is
     a database focused on radio astronomy data and images, primarily
     those collected by ASTRON itself.  A query to ‘astron’ is submitted
     to ‘https://vo.astron.nl/__system__/tap/run/tap/sync’.

     Here is the list of short names for dataset(s) in ASTRON's VO
     service:
        • ‘tgssadr --> tgssadr.main’

‘gaia’
     The Gaia project (<https://www.cosmos.esa.int/web/gaia>) database
     which is a large collection of star positions on the celestial
     sphere, as well as peculiar velocities, parallaxes and magnitudes
     in some bands among many others.  Besides scientific studies (like
     studying resolved stellar populations in the Galaxy and its halo),
     Gaia is also invaluable for raw data calibrations, like astrometry.
     A query to ‘gaia’ is submitted to
     ‘https://gea.esac.esa.int/tap-server/tap/sync’.

     Here is the list of short names for popular datasets within Gaia:
        • ‘dr3 --> gaiadr3.gaia_source’
        • ‘edr3 --> gaiaedr3.gaia_source’
        • ‘dr2 --> gaiadr2.gaia_source’
        • ‘dr1 --> gaiadr1.gaia_source’
        • ‘tycho2 --> public.tycho2’
        • ‘hipparcos --> public.hipparcos’

‘ned’
     The NASA/IPAC Extragalactic Database (NED,
     <http://ned.ipac.caltech.edu>) is a fusion database, integrating
     the information about extra-galactic sources from many large sky
     surveys into a single catalog.  It covers the full spectrum, from
     Gamma rays to radio frequencies and is updated when new data
     arrives.  A TAP query to ‘ned’ is submitted to
     ‘https://ned.ipac.caltech.edu/tap/sync’.

        • ‘objdir --> NEDTAP.objdir’: default TAP-based dataset in NED.

        • ‘extinction’: A command-line interface to the NED Extinction
          Calculator
          (https://ned.ipac.caltech.edu/extinction_calculator).  It only
          takes a central coordinate and returns a VOTable of the
          calculated extinction in many commonly used filters at that
          point.  As a result, options like ‘--width’ or ‘--radius’ are
          not supported.  However, Gnuastro does not yet support the
          VOTable format.  Therefore, if you specify an ‘--output’ file,
          it should have an ‘.xml’ suffix and the downloaded file will
          not be checked.

          Until VOTable support is added to Gnuastro, you can use GREP,
          AWK and SED to convert the VOTable data into a FITS table with
          a command like below (assuming the queried VOTable is called
          ‘ned-extinction.xml’):

          grep '^<TR><TD>' ned-extinction.xml \
           | sed -e's|<TR><TD>||' \
                 -e's|</TD></TR>||' \
                 -e's|</TD><TD>|@|g' \
           | awk 'BEGIN{FS="@"; \
               print "# Column 1: FILTER [name,str15] Filter name"; \
               print "# Column 2: CENTRAL [um,f32] Central Wavelength"; \
               print "# Column 3: EXTINCTION [mag,f32] Galactic Ext."; \
               print "# Column 4: ADS_REF [ref,str50] ADS reference"} \
                  {printf "%-15s %g %g %s\n", $1, $2, $3, $4}' \
           | asttable -oned-extinction.fits

          Once the table is in FITS, you can easily get the extinction
          for a certain filter (for example, the ‘SDSS r’ filter) like
          the command below:

               asttable ned-extinction.fits --equal=FILTER,"SDSS r" \
                        -cEXTINCTION

‘vizier’
     Vizier (<https://vizier.u-strasbg.fr>) is arguably the largest
     catalog database in astronomy: containing more than 20500 catalogs
     as of mid January 2021.  Almost all published catalogs in major
     projects, and even the tables in many papers are archived and
     accessible here.  For example, VizieR also has a full copy of the
     Gaia database mentioned below, with some additional standardized
     columns (like RA and Dec in J2000).

     The current implementation of ‘--limitinfo’ only looks into the
     description of the datasets, but since VizieR is so large, there is
     still a lot of room for improvement.  Until then, if ‘--limitinfo’
     is not sufficient, you can use VizieR's own web-based search for
     your desired dataset: <http://cdsarc.u-strasbg.fr/viz-bin/cat>

     Because VizieR curates such a diverse set of data from tens of
     thousands of projects and aims for interoperability between them,
     the column names in VizieR may not be identical to the column names
     in the surveys' own databases (Gaia in the example above).  A query
     to ‘vizier’ is submitted to
     ‘http://tapvizier.u-strasbg.fr/TAPVizieR/tap/sync’.

     Here is the list of short names for popular datasets within VizieR
     (sorted alphabetically by their short name).  Please feel free to
     suggest other major catalogs (covering a wide area or commonly used
     in your field)..  For details on each dataset with necessary
     citations, and links to web pages, look into their details with
     their ViziR names in <https://vizier.u-strasbg.fr/viz-bin/VizieR>.
        • ‘2mass --> II/246/out’ (2MASS All-Sky Catalog)
        • ‘akarifis --> II/298/fis’ (AKARI/FIS All-Sky Survey)
        • ‘allwise --> II/328/allwise’ (AllWISE Data Release)
        • ‘apass9 --> II/336/apass9’ (AAVSO Photometric All Sky Survey,
          DR9)
        • ‘catwise --> II/365/catwise’ (CatWISE 2020 catalog)
        • ‘des1 --> II/357/des_dr1’ (Dark Energy Survey data release 1)
        • ‘gaiadr3 --> I/355/gaiadr3’ (GAIA Data Release 3)
        • ‘gaiaedr3 --> I/350/gaiadr3’ (GAIA Early Data Release 3)
        • ‘gaiadr2 --> I/345/gaia2’ (GAIA Data Release 2)
        • ‘galex5 --> II/312/ais’ (All-sky Survey of GALEX DR5)
        • ‘nomad --> I/297/out’ (Naval Observatory Merged Astrometric
          Dataset)
        • ‘panstarrs1 --> II/349/ps1’ (Pan-STARRS Data Release 1).
        • ‘ppmxl --> I/317/sample’ (Positions and proper motions on the
          ICRS)
        • ‘sdss12 --> V/147/sdss12’ (SDSS Photometric Catalogue, Release
          12)
        • ‘usnob1 --> I/284/out’ (Whole-Sky USNO-B1.0 Catalog)
        • ‘ucac5 --> I/340/ucac5’ (5th U.S. Naval Obs.  CCD Astrograph
          Catalog)
        • ‘unwise --> II/363/unwise’ (Band-merged unWISE Catalog)
        • ‘wise --> II/311/wise’ (WISE All-Sky data Release)


File: gnuastro.info,  Node: Invoking astquery,  Prev: Available databases,  Up: Query

5.4.2 Invoking Query
--------------------

Query provides a high-level interface to downloading subsets of data
from databases.  The executable name is ‘astquery’ with the following
general template

     $ astquery DATABASE-NAME [OPTION...] ...

One line examples:


     ## Information about all datasets in ESA's GAIA database:
     $ astquery gaia --information

     ## Only show catalogs in VizieR that have 'MUSE' in their
     ## description. The '-i' is short for '--information'.
     $ astquery vizier -i --limitinfo=MUSE

     ## List of columns in 'J/A+A/608/A2/udf10' (one of the above).
     $ astquery vizier --dataset=J/A+A/608/A2/udf10 -i

     ## ID, RA and Dec of all Gaia sources within an image.
     $ astquery gaia --dataset=dr3 --overlapwith=image.fits \
                -csource_id,ra,dec

     ## RA, Dec and Spectroscopic redshifts of objects in SDSS DR12
     ## spectroscopic redshift that overlap with 'image.fits'.
     $ astquery vizier --dataset=sdss12 --overlapwith=image.fits \
                -cRA_ICRS,DE_ICRS,zsp --range=zsp,1e-10,inf

     ## All columns of all entries in the Gaia DR3 catalog (hosted at
     ## VizieR) within 1 arc-minute of the given coordinate.
     $ astquery vizier --dataset=gaiadr3 --output=my-gaia.fits \
                --center=113.8729761,31.9027152 --radius=1/60 \

     ## Similar to above, but only ID, RA and Dec columns for objects with
     ## magnitude range 10 to 15. In VizieR, this column is called 'Gmag'.
     ## Also, using sexagesimal coordinates instead of degrees for center.
     $ astquery vizier --dataset=gaiadr3 --output=my-gaia.fits \
                --center=07h35m29.51,31d54m9.77 --radius=1/60 \
                --range=Gmag,10:15 -cDR3Name,RAJ2000,DEJ2000

   Query takes a single argument which is the name of the database.  For
the full list of available databases and accessing them, see *note
Available databases::.  There are two methods to query the databases,
each is more fully discussed in its option's description below.
   • *Low-level:* With ‘--query’ you can directly give a raw query
     statement that is recognized by the database.  This is very low
     level and will require a good knowledge of the database's query
     language, but of course, it is much more powerful.  If this option
     is given, the raw string is directly passed to the server and all
     other constraints/options (for Query's high-level interface) are
     ignored.
   • *High-level:* With the high-level options (like ‘--column’,
     ‘--center’, ‘--radius’, ‘--range’ and other constraining options
     below), the low-level query will be constructed automatically for
     the particular database.  This method is only limited to the
     generic capabilities that Query provides for all servers.  So
     ‘--query’ is more powerful, however, in this mode, you do not need
     any knowledge of the database's query language.  You can see the
     internally generated query on the terminal (if ‘--quiet’ is not
     used) or in the 0-th extension of the output (if it is a FITS
     file).  This full command contains the internally generated query.

   The name of the downloaded output file can be set with ‘--output’.
The requested output format can have any of the *note Recognized table
formats:: (currently ‘.txt’ or ‘.fits’).  Like all Gnuastro programs, if
the output is a FITS file, the zero-th/first HDU of the output will
contain all the command-line options given to Query as well as the full
command used to access the server.  When ‘--output’ is not set, the
output name will be in the format of ‘NAME-STRING.fits’, where ‘NAME’ is
the name of the database and ‘STRING’ is a randomly selected 6-character
set of numbers and alphabetic characters.  With this feature, a second
run of ‘astquery’ that is not called with ‘--output’ will not over-write
an already downloaded one.  Generally, when calling Query more than
once, it is recommended to set an output name for each call based on
your project's context.

   The outputs of Query will have a common output format, irrespective
of the used database.  To achieve this, Query will ask the databases to
provide a FITS table output (for larger tables, FITS can consume much
less download volume).  After downloading is complete, the raw
downloaded file will be read into memory once by Query, and written into
the file given to ‘--output’.  The raw downloaded file will be deleted
by default, but can be preserved with the ‘--keeprawdownload’ option.
This strategy avoids unnecessary surprises depending on database.  For
example, some databases can download a compressed FITS table, even
though we ask for FITS. But with the strategy above, the final output
will be an uncompressed FITS file.  The metadata that is added by Query
(including the full download command) is also very useful for future
usage of the downloaded data.  Unfortunately many databases do not write
the input queries into their generated tables.

‘--dry-run’
     Only print the final download command to contact the server, do not
     actually run it.  This option is good when you want to check the
     finally constructed query or download options given to the download
     program.  You may also want to use the constructed command as a
     base to do further customizations on it and run it yourself.

‘-k’
‘--keeprawdownload’
     Do not delete the raw downloaded file from the database.  The name
     of the raw download will have a ‘OUTPUT-raw-download.fits’ format.
     Where ‘OUTPUT’ is either the base-name of the final output file
     (without a suffix).

‘-i’
‘--information’
     Print the information of all datasets (tables) within a database or
     all columns within a database.  When ‘--dataset’ is specified, the
     latter mode (all column information) is downloaded and printed and
     when it is not defined, all dataset information (within the
     database) is printed.

     Some databases (like VizieR) contain tens of thousands of datasets,
     so you can limit the downloaded and printed information for
     available databases with the ‘--limitinfo’ option (described
     below).  Dataset descriptions are often large and contain a lot of
     text (unlike column descriptions).  Therefore when printing the
     information of all datasets within a database, the information
     (e.g., database name) will be printed on separate lines before the
     description.  However, when printing column information, the output
     has the same format as a similar option in Table (see *note
     Invoking asttable::).

     Important note to consider: the printed order of the datasets or
     columns is just for displaying in the printed output.  You cannot
     ask for datasets or columns based on the printed order, you need to
     use dataset or column names.

‘-L STR’
‘--limitinfo=STR’
     Limit the information that is downloaded and displayed (with
     ‘--information’) to those that have the string given to this option
     in their description.  Note that _this is case-sensitive_.  This
     option is only relevant when ‘--information’ is also called.

     Databases may have thousands (or tens of thousands) of datasets.
     Therefore just the metadata (information) to show with
     ‘--information’ can be tens of megabytes (for example, the full
     VizieR metadata file is about 23Mb as of January 2021).  Once
     downloaded, it can also be hard to parse manually.  With
     ‘--limitinfo’, only the metadata of datasets that contain this
     string _in their description_ will be downloaded and displayed,
     greatly improving the speed of finding your desired dataset.

‘-Q "STR"’
‘--query="STR"’
     Directly specify the query to be passed onto the database.  The
     queries will generally contain space and other meta-characters, so
     we recommend placing the query within quotations.

‘-s STR’
‘--dataset=STR’
     The dataset to query within the database (not compatible with
     ‘--query’).  This option is mandatory when ‘--query’ or
     ‘--information’ are not provided.  You can see the list of
     available datasets within a database using ‘--information’
     (possibly supplemented by ‘--limitinfo’).  The output of
     ‘--information’ will contain the recognized name of the datasets
     within that database.  You can pass the recognized name directly to
     this option.  For more on finding and using your desired database,
     see *note Available databases::.

‘-c STR’
‘--column=STR[,STR[,...]]’
     The column name(s) to retrieve from the dataset in the given order
     (not compatible with ‘--query’).  If not given, all the dataset's
     columns for the selected rows will be queried (which can be
     large!).  This option can take multiple values in one instance (for
     example, ‘--column=ra,dec,mag’), or in multiple instances (for
     example, ‘-cra -cdec -cmag’), or mixed (for example, ‘-cra,dec
     -cmag’).

     In case, you do not know the full list of the dataset's column
     names a-priori, and you do not want to download all the columns
     (which can greatly decrease your download speed), you can use the
     ‘--information’ option combined with the ‘--dataset’ option, see
     *note Available databases::.

‘-H INT’
‘--head=INT’
     Only ask for the first ‘INT’ rows of the finally selected columns,
     not all the rows.  This can be good when your search can result a
     large dataset, but before downloading the full volume, you want to
     see the top rows and get a feeling of what the whole dataset looks
     like.

‘-v FITS’
‘--overlapwith=FITS’
     File name of FITS file containing an image (in the HDU given by
     ‘--hdu’) to use for identifying the region to query in the give
     database and dataset.  Based on the image's WCS and pixel size, the
     sky coverage of the image is estimated and values to the
     ‘--center’, ‘--width’ will be calculated internally.  Hence this
     option cannot be used with ‘--center’, ‘--width’ or ‘--radius’.
     Also, since it internally generates the query, it cannot be used
     with ‘--query’.

     If the image is rotated in relation to RA/DEC, or the image has a
     large coverage on the sky, or it has WCS distortions (that the
     server can't know about), the retrieved catalog may not fully
     overlap with the image (correspond to a larger area in the sky).
     To be sure that the final catalog you use actually has sources
     within the image, use the commands below.  Two points to have in
     mind in the example below: 1) use the ‘cols’ variable to specify
     the names of the columns you want (this is necessary since this
     list of column names is necessary in two places of the commands
     below).  2) The last two commands are just for visual validation,
     they are not necessary in a script.

          $ cols=source_id,ra,dec,phot_g_mean_mag

          $ astquery gaia --dataset=dr3 --overlapwith=image.fits \
                     --column=$cols --output=gaia.fits

          $ astcrop image.fits --catalog=gaia.fits --coordcol=ra \
                    --coordcol=dec --mode=wcs --width=1 \
                    --widthinpix --log --oneelemstdout --quiet \
                    --output=crop-log.fits > /dev/null

          $ asttable gaia.fits --catcolumnfile=crop-log.fits \
                     --catcolumns=NUM_INPUTS -O \
                     | asttable --equal=NUM_INPUTS,1 \
                                --column=$cols \
                                --output=gaia-in-image.fits

          $ astscript-ds9-region gaia-in-image.fits --column=ra,dec \
                                 --radius=5 --width=5 \
                                 --output=stars-in.reg

          $ astscript-fits-view image.fits --region=stars-in.reg

     Note that if the image has WCS distortions and the reference point
     for the WCS is not within the image, the WCS will not be
     well-defined.  Therefore the resulting catalog may not overlap, or
     correspond to a larger/small area in the sky.

‘-C FLT,FLT’
‘--center=FLT,FLT’
     The spatial center position (mostly RA and Dec) to use for the
     automatically generated query (not compatible with ‘--query’).  The
     comma-separated values can either be in degrees (a single number),
     or sexagesimal (‘_h_m_’ for RA, ‘_d_m_’ for Dec, or ‘_:_:_’ for
     both).

     The given values will be compared to two columns in the database to
     find/return rows within a certain region around this center
     position will be requested and downloaded.  Pre-defined RA and Dec
     column names are defined in Query for every database, however you
     can use ‘--ccol’ to select other columns to use instead.  The
     region can either be a circle and the point (configured with
     ‘--radius’) or a box/rectangle around the point (configured with
     ‘--width’).

‘--ccol=STR,STR’
     The name of the coordinate-columns in the dataset to compare with
     the values given to ‘--center’.  Query will use its internal
     defaults for each dataset (for example, ‘RAJ2000’ and ‘DEJ2000’ for
     VizieR data).  But each dataset is treated separately and it is not
     guaranteed that these columns exist in all datasets.  Also, more
     than one coordinate system/epoch may be present in a dataset and
     you can use this option to construct your spatial constraint based
     on the others coordinate systems/epochs.

‘-r FLT’
‘--radius=FLT’
     The radius about the requested center to use for the automatically
     generated query (not compatible with ‘--query’).  The radius is in
     units of degrees, but you can use simple division with this option
     directly on the command-line.  For example, if you want a radius of
     20 arc-minutes or 20 arc-seconds, you can use ‘--radius=20/60’ or
     ‘--radius=20/3600’ respectively (which is much more human-friendly
     than ‘0.3333’ or ‘0.005556’).

‘-w FLT[,FLT]’
‘--width=FLT[,FLT]’
     The square (or rectangle) side length (width) about the requested
     center to use for the automatically generated query (not compatible
     with ‘--query’).  If only one value is given to ‘--width’ the
     region will be a square, but if two values are given, the widths of
     the query box along each dimension will be different.  The value(s)
     is (are) in the same units as the coordinate column (see ‘--ccol’,
     usually RA and Dec which are degrees).  You can use simple division
     for each value directly on the command-line if you want relatively
     small (and more human-friendly) sizes.  For example, if you want
     your box to be 1 arc-minutes along the RA and 2 arc-minutes along
     Dec, you can use ‘--width=1/60,2/60’.

‘-g STR,FLT,FLT’
‘--range=STR,FLT,FLT’
     The column name and numerical range (inclusive) of acceptable
     values in that column (not compatible with ‘--query’).  This option
     can be called multiple times for applying range limits on many
     columns in one call (thus greatly reducing the download size).  For
     example, when used on the ESA gaia database, you can use
     ‘--range=phot_g_mean_mag,10:15’ to only get rows that have a value
     between 10 and 15 (inclusive on both sides) in the
     ‘phot_g_mean_mag’ column.

     If you want all rows larger, or smaller, than a certain number, you
     can use ‘inf’, or ‘-inf’ as the first or second values
     respectively.  For example, if you want objects with SDSS
     spectroscopic redshifts larger than 2 (from the VizieR ‘sdss12’
     database), you can use ‘--range=zsp,2,inf’

     If you want the interval to not be inclusive on both sides, you can
     run ‘astquery’ once and get the command that it executes.  Then you
     can edit it to be non-inclusive on your desired side.

‘-b STR[,STR]’
‘--noblank=STR[,STR]’
     Only ask for rows that do not have a blank value in the ‘STR’
     column.  This option can be called many times, and each call can
     have multiple column names (separated by a comma or <,>).  For
     example, if you want the retrieved rows to not have a blank value
     in columns ‘A’, ‘B’, ‘C’ and ‘D’, you can use ‘--noblank=A
     -bB,C,D’.

‘--sort=STR[,STR]’
     Ask for the server to sort the downloaded data based on the given
     columns.  For example, let's assume your desired catalog has column
     ‘Z’ for redshift and column ‘MAG_R’ for magnitude in the R band.
     When you call ‘--sort=Z,MAG_R’, it will primarily sort the columns
     based on the redshift, but if two objects have the same redshift,
     they will be sorted by magnitude.  You can add as many columns as
     you like for higher-level sorting.


File: gnuastro.info,  Node: Data manipulation,  Next: Data analysis,  Prev: Data containers,  Up: Top

6 Data manipulation
*******************

Images are one of the major formats of data that is used in astronomy.
The functions in this chapter explain the GNU Astronomy Utilities which
are provided for their manipulation.  For example, cropping out a part
of a larger image or convolving the image with a given kernel or
applying a transformation to it.

* Menu:

* Crop::                        Crop region(s) from a dataset.
* Arithmetic::                  Arithmetic on input data.
* Convolve::                    Convolve an image with a kernel.
* Warp::                        Warp/Transform an image to a different grid.


File: gnuastro.info,  Node: Crop,  Next: Arithmetic,  Prev: Data manipulation,  Up: Data manipulation

6.1 Crop
========

Astronomical images are often very large, filled with thousands of
galaxies.  It often happens that you only want a section of the image,
or you have a catalog of sources and you want to visually analyze them
in small postage stamps.  Crop is made to do all these things.  When
more than one crop is required, Crop will divide the crops between
multiple threads to significantly reduce the run time.

   Astronomical surveys are usually extremely large.  So large in fact,
that the whole survey will not fit into a reasonably sized file.
Because of this, surveys usually cut the final image into separate tiles
and store each tile in a file.  For example, the COSMOS survey's Hubble
space telescope, ACS F814W image consists of 81 separate FITS images,
with each one having a volume of 1.7 Gigabytes.

   Even though the tile sizes are chosen to be large enough that too
many galaxies/targets do not fall on the edges of the tiles, inevitably
some do.  So when you simply crop the image of such targets from one
tile, you will miss a large area of the surrounding sky (which is
essential in estimating the noise).  Therefore in its WCS mode, Crop
will stitch parts of the tiles that are relevant for a target (with the
given width) from all the input images that cover that region into the
output.  Of course, the tiles have to be present in the list of input
files.

   Besides cropping postage stamps around certain coordinates, Crop can
also crop arbitrary polygons from an image (or a set of tiles by
stitching the relevant parts of different tiles within the polygon), see
‘--polygon’ in *note Invoking astcrop::.  Alternatively, it can crop out
rectangular regions through the ‘--section’ option from one image, see
*note Crop section syntax::.

* Menu:

* Crop modes::                  Basic modes to define crop region.
* Crop section syntax::         How to define a section to crop.
* Blank pixels::                Pixels with no value.
* Invoking astcrop::            Calling Crop on the command-line


File: gnuastro.info,  Node: Crop modes,  Next: Crop section syntax,  Prev: Crop,  Up: Crop

6.1.1 Crop modes
----------------

In order to be comprehensive, intuitive, and easy to use, there are two
ways to define the crop:

  1. From its center and side length.  For example, if you already know
     the coordinates of an object and want to inspect it in an image or
     to generate postage stamps of a catalog containing many such
     coordinates.

  2. The vertices of the crop region, this can be useful for larger
     crops over many targets, for example, to crop out a uniformly deep,
     or contiguous, region of a large survey.

   Irrespective of how the crop region is defined, the coordinates to
define the crop can be in Image (pixel) or World Coordinate System (WCS)
standards.  All coordinates are read as floating point numbers (not
integers, except for the ‘--section’ option, see below).  By setting the
_mode_ in Crop, you define the standard that the given coordinates must
be interpreted.  Here, the different ways to specify the crop region are
discussed within each standard.  For the full list options, please see
*note Invoking astcrop::.

   When the crop is defined by its center, the respective (integer)
central pixel position will be found internally according to the FITS
standard.  To have this pixel positioned in the center of the cropped
region, the final cropped region will have an add number of pixels (even
if you give an even number to ‘--width’ in image mode).

   Furthermore, when the crop is defined as by its center, Crop allows
you to only keep crops what do not have any blank pixels in the vicinity
of their center (your primary target).  This can be very convenient when
your input catalog/coordinates originated from another survey/filter
which is not fully covered by your input image, to learn more about this
feature, please see the description of the ‘--checkcenter’ option in
*note Invoking astcrop::.

Image coordinates
     In image mode (‘--mode=img’), Crop interprets the pixel coordinates
     and widths in units of the input data-elements (for example, pixels
     in an image, not world coordinates).  In image mode, only one image
     may be input.  The output crop(s) can be defined in multiple ways
     as listed below.

     Center of multiple crops (in a catalog)
          The center of (possibly multiple) crops are read from a text
          file.  In this mode, the columns identified with the
          ‘--coordcol’ option are interpreted as the center of a crop
          with a width of ‘--width’ pixels along each dimension.  The
          columns can contain any floating point value.  The value to
          ‘--output’ option is seen as a directory which will host (the
          possibly multiple) separate crop files, see *note Crop
          output:: for more.  For a tutorial using this feature, please
          see *note Reddest clumps cutouts and parallelization::.

     Center of a single crop (on the command-line)
          The center of the crop is given on the command-line with the
          ‘--center’ option.  The crop width is specified by the
          ‘--width’ option along each dimension.  The given coordinates
          and width can be any floating point number.

     Vertices of a single crop
          In Image mode there are two options to define the vertices of
          a region to crop: ‘--section’ and ‘--polygon’.  The former is
          lower-level (does not accept floating point vertices, and only
          a rectangular region can be defined), it is also only
          available in Image mode.  Please see *note Crop section
          syntax:: for a full description of this method.

          The latter option (‘--polygon’) is a higher-level method to
          define any polygon (with any number of vertices) with floating
          point values.  Please see the description of this option in
          *note Invoking astcrop:: for its syntax.

WCS coordinates
     In WCS mode (‘--mode=wcs’), the coordinates and width are
     interpreted using the World Coordinate System (WCS, that must
     accompany the dataset), not pixel coordinates.  You can optionally
     use ‘--widthinpix’ for the width to be interpreted in pixels (even
     though the coordinates are in WCS). In WCS mode, Crop accepts
     multiple datasets as input.  When the cropped region (defined by
     its center or vertices) overlaps with multiple of the input
     images/tiles, the overlapping regions will be taken from the
     respective input (they will be stitched when necessary for each
     output crop).

     In this mode, the input images do not necessarily have to be the
     same size, they just need to have the same orientation and pixel
     resolution.  Currently only orientation along the celestial
     coordinates is accepted, if your input has a different orientation
     or resolution you can use Warp's ‘--gridfile’ option to align the
     image before cropping it (see *note Warp::).

     Each individual input image/tile can even be smaller than the final
     crop.  In any case, any part of any of the input images which
     overlaps with the desired region will be used in the crop.  Note
     that if there is an overlap in the input images/tiles, the pixels
     from the last input image read are going to be used for the
     overlap.  Crop will not change pixel values, so it assumes your
     overlapping tiles were cutout from the same original image.  There
     are multiple ways to define your cropped region as listed below.

     Center of multiple crops (in a catalog)
          Similar to catalog inputs in Image mode (above), except that
          the values along each dimension are assumed to have the same
          units as the dataset's WCS information.  For example, the
          central RA and Dec value for each crop will be read from the
          first and second calls to the ‘--coordcol’ option.  The width
          of the cropped box (in units of the WCS, or degrees in RA and
          Dec mode) must be specified with the ‘--width’ option.  You
          can optionally use ‘--widthinpix’ for the value of ‘--width’
          to be interpreted in pixels.

     Center of a single crop (on the command-line)
          You can specify the center of only one crop box with the
          ‘--center’ option.  If it exists in the input images, it will
          be cropped similar to the catalog mode, see above also for
          ‘--width’.

     Vertices of a single crop
          The ‘--polygon’ option is a high-level method to define any
          convex polygon (with any number of vertices).  Please see the
          description of this option in *note Invoking astcrop:: for its
          syntax.

     *CAUTION:* In WCS mode, the image has to be aligned with the
     celestial coordinates, such that the first FITS axis is parallel
     (opposite direction) to the Right Ascension (RA) and the second
     FITS axis is parallel to the declination.  If these conditions are
     not met for an image, Crop will warn you and abort.  You can use
     Warp to align the input image to standard celestial coordinates,
     see *note Warp::.

   As a summary, if you do not specify a catalog, you have to define the
cropped region manually on the command-line.  In any case the mode is
mandatory for Crop to be able to interpret the values given as
coordinates or widths.


File: gnuastro.info,  Node: Crop section syntax,  Next: Blank pixels,  Prev: Crop modes,  Up: Crop

6.1.2 Crop section syntax
-------------------------

When in image mode, one of the methods to crop only one rectangular
section from the input image is to use the ‘--section’ option.  Crop has
a powerful syntax to read the box parameters from a string of
characters.  If you leave certain parts of the string to be empty, Crop
can fill them for you based on the input image sizes.

   To define a box, you need the coordinates of two points: the first
(‘X1’, ‘Y1’) and the last pixel (‘X2’, ‘Y2’) pixel positions in the
image, or four integer numbers in total.  The four coordinates can be
specified with one string in this format: '‘X1:X2,Y1:Y2’'.  This string
is given to the ‘--section’ option.  Therefore, the pixels along the
first axis that are $\geq$‘X1’ and $\leq$‘X2’ will be included in the
cropped image.  The same goes for the second axis.  Note that each
different term will be read as an integer, not a float.

   The reason it only accepts integers is that ‘--section’ is a
low-level option (which is also very fast!).  For a higher-level way to
specify region (any polygon, not just a box), please see the ‘--polygon’
option in *note Crop options::.  Also note that in the FITS standard,
pixel indexes along each axis start from unity(1) not zero(0).

   You can omit any of the values and they will be filled automatically.
The left hand side of the colon (‘:’) will be filled with ‘1’, and the
right side with the image size.  So, ‘2:,:’ will include the full range
of pixels along the second axis and only those with a first axis index
larger than ‘2’ in the first axis.  If the colon is omitted for a
dimension, then the full range is automatically used.  So the same
string is also equal to ‘2:,’ or ‘2:’ or even ‘2’.  If you want such a
case for the second axis, you should set it to: ‘,2’.

   If you specify a negative value, it will be seen as before the
indexes of the image which are outside the image along the bottom or
left sides when viewed in SAO DS9.  In case you want to count from the
top or right sides of the image, you can use an asterisk (‘*’).  When
confronted with a ‘*’, Crop will replace it with the maximum length of
the image in that dimension.  So ‘*-10:*+10,*-20:*+20’ will mean that
the crop box will be 20\times40 pixels in size and only include the top
corner of the input image with 3/4 of the image being covered by blank
pixels, see *note Blank pixels::.

   If you feel more comfortable with space characters between the
values, you can use as many space characters as you wish, just be
careful to put your value in double quotes, for example,
‘--section="5:200, 123:854"’.  If you forget the quotes, anything after
the first space will not be seen by ‘--section’ and you will most
probably get an error because the rest of your string will be read as a
filename (which most probably does not exist).  See *note Command-line::
for a description of how the command-line works.


File: gnuastro.info,  Node: Blank pixels,  Next: Invoking astcrop,  Prev: Crop section syntax,  Up: Crop

6.1.3 Blank pixels
------------------

The cropped box can potentially include pixels that are beyond the image
range.  For example, when a target in the input catalog was very near
the edge of the input image.  The parts of the cropped image that were
not in the input image will be filled with the following two values
depending on the data type of the image.  In both cases, SAO DS9 will
not color code those pixels.
   • If the data type of the image is a floating point type (float or
     double), IEEE NaN (Not a number) will be used.
   • For integer types, pixels out of the image will be filled with the
     value of the ‘BLANK’ keyword in the cropped image header.  The
     value assigned to it is the lowest value possible for that type, so
     you will probably never need it any way.  Only for the unsigned
     character type (‘BITPIX=8’ in the FITS header), the maximum value
     is used because it is unsigned, the smallest value is zero which is
     often meaningful.
   You can ask for such blank regions to not be included in the output
crop image using the ‘--noblank’ option.  In such cases, there is no
guarantee that the image size of your outputs are what you asked for.

   In some survey images, unfortunately they do not use the ‘BLANK’ FITS
keyword.  Instead they just give all pixels outside of the survey area a
value of zero.  So by default, when dealing with float or double image
types, any values that are 0.0 are also regarded as blank regions.  This
can be turned off with the ‘--zeroisnotblank’ option.


File: gnuastro.info,  Node: Invoking astcrop,  Prev: Blank pixels,  Up: Crop

6.1.4 Invoking Crop
-------------------

Crop will crop a region from an image.  If in WCS mode, it will also
stitch parts from separate images in the input files.  The executable
name is ‘astcrop’ with the following general template

     $ astcrop [OPTION...] [ASCIIcatalog] ASTRdata ...

One line examples:

     ## Crop all objects in cat.txt from image.fits:
     $ astcrop --catalog=cat.txt image.fits

     ## Crop all options in catalog (with RA,DEC) from all the files
     ## ending in `_drz.fits' in `/mnt/data/COSMOS/':
     $ astcrop --mode=wcs --catalog=cat.txt /mnt/data/COSMOS/*_drz.fits

     ## Crop the outer 10 border pixels of the input image and give
     ## the output HDU a name ('EXTNAME' keyword in FITS) of 'mysection'.
     $ astcrop --section=10:*-10,10:*-10 --hdu=2 image.fits \
               --metaname=mysection

     ## Crop region around RA and Dec of (189.16704, 62.218203):
     $ astcrop --mode=wcs --center=189.16704,62.218203 goodsnorth.fits

     ## Same crop above, but coordinates given in sexagesimal (you can
     ## also use ':' between the sexagesimal components).
     $ astcrop --mode=wcs --center=12h36m40.08,62d13m5.53 goodsnorth.fits

     ## Crop region around pixel coordinate (568.342, 2091.719):
     $ astcrop --mode=img --center=568.342,2091.719 --width=201 image.fits

     ## Crop all HDUs within a FITS file at a certain coordinate, while
     ## preserving the names of the HDUs in the output.
     $ for hdu in $(astfits input.fits --listimagehdus); do \
         astcrop input.fits --hdu=$hdu --append --output=crop.fits \
                 --metaname=$hdu --mode=wcs --center=189.16704,62.218203 \
                 --width=10/3600
       done

Crop has one mandatory argument which is the input image name(s), shown
above with ‘ASTRdata ...’.  You can use shell expansions, for example,
‘*’ for this if you have lots of images in WCS mode.  If the crop box
centers are in a catalog, you can use the ‘--catalog’ option.  In other
cases, you have to provide the single cropped output parameters must be
given with command-line options.  See *note Crop output:: for how the
output file name(s) can be specified.  For the full list of general
options to all Gnuastro programs (including Crop), please see *note
Common options::.

   Floating point numbers can be used to specify the crop region (except
the ‘--section’ option, see *note Crop section syntax::).  In such
cases, the floating point values will be used to find the desired
integer pixel indices based on the FITS standard.  Hence, Crop
ultimately does not do any sub-pixel cropping (in other words, it does
not change pixel values).  If you need such crops, you can use *note
Warp:: to first warp the image to the a new pixel grid, then crop from
that.  For example, let's assume you want a crop from pixels 12.982 to
80.982 along the first dimension.  You should first translate the image
by $-0.482$ (note that the edge of a pixel is at integer multiples of
$0.5$).  So you should run Warp with ‘--translate=-0.482,0’ and then
crop the warped image with ‘--section=13:81’.

   There are two ways to define the cropped region: with its center or
its vertices.  See *note Crop modes:: for a full description.  In the
former case, Crop can check if the central region of the cropped image
is indeed filled with data or is blank (see *note Blank pixels::), and
not produce any output when the center is blank, see the description
under ‘--checkcenter’ for more.

   When in catalog mode, Crop will run in parallel unless you set
‘--numthreads=1’, see *note Multi-threaded operations::.  Note that when
multiple outputs are created with threads, the outputs will not be
created in the same order.  This is because the threads are asynchronous
and thus not started in order.  This has no effect on each output, see
*note Reddest clumps cutouts and parallelization:: for a tutorial on
effectively using this feature.

* Menu:

* Crop options::                A list of all the options with explanation.
* Crop output::                 The outputs of Crop.
* Crop known issues::           Known issues in running Crop.


File: gnuastro.info,  Node: Crop options,  Next: Crop output,  Prev: Invoking astcrop,  Up: Invoking astcrop

6.1.4.1 Crop options
....................

The options can be classified into the following contexts: Input, Output
and operating mode options.  Options that are common to all Gnuastro
program are listed in *note Common options:: and will not be repeated
here.

   When you are specifying the crop vertices yourself (through
‘--section’, or ‘--polygon’) on relatively small regions (depending on
the resolution of your images) the outputs from image and WCS mode can
be approximately equivalent.  However, as the crop sizes get large, the
curved nature of the WCS coordinates have to be considered.  For
example, when using ‘--section’, the right ascension of the bottom left
and top left corners will not be equal.  If you only want regions within
a given right ascension, use ‘--polygon’ in WCS mode.

Input image parameters:

‘--hstartwcs=INT’
     Specify the first keyword card (line number) to start finding the
     input image world coordinate system information.  This is useful
     when certain header keywords of the input may cause bad conflicts
     with your crop (see an example described below).  To get line
     numbers of the header keywords, you can pipe the fully printed
     header into ‘cat -n’ like below:

          $ astfits image.fits -h1 | cat -n

     For example, distortions have only been present in WCSLIB from
     version 5.15 (released in mid 2016).  Therefore some pipelines
     still apply their own specific set of WCS keywords for distortions
     and put them into the image header along with those that WCSLIB
     does recognize.  So now that WCSLIB recognizes most of the standard
     distortion parameters, they will get confused with the old ones and
     give wrong results.  For example, in the CANDELS-GOODS South images
     that were created before WCSLIB 5.15(1).

     The two ‘--hstartwcs’ and ‘--hendwcs’ are thus provided so when
     using older datasets, you can specify what region in the FITS
     headers you want to use to read the WCS keywords.  Note that this
     is only relevant for reading the WCS information, basic data
     information like the image size are read separately.  These two
     options will only be considered when the value to ‘--hendwcs’ is
     larger than that of ‘--hstartwcs’.  So if they are equal or
     ‘--hstartwcs’ is larger than ‘--hendwcs’, then all the input
     keywords will be parsed to get the WCS information of the image.

‘--hendwcs=INT’
     Specify the last keyword card to read for specifying the image
     world coordinate system on the input images.  See ‘--hstartwcs’

Crop box parameters:

‘-c FLT[,FLT[,...]]’
‘--center=FLT[,FLT[,...]]’
     The central position of the crop in the input image.  The positions
     along each dimension must be separated by a comma (<,>) and
     fractions are also acceptable.  The comma-separated values can
     either be in degrees (a single number), or sexagesimal (‘_h_m_’ for
     RA, ‘_d_m_’ for Dec, or ‘_:_:_’ for both).

     The number of values given to this option must be the same as the
     dimensions of the input dataset.  The width of the crop should be
     set with ‘--width’.  The units of the coordinates are read based on
     the value to the ‘--mode’ option, see below.

‘-O STR’
‘--mode=STR’
     Mode to interpret the crop's coordinates (for example with
     ‘--center’, ‘--catalog’ or ‘--polygon’).  The value must either be
     ‘img’ (to assume image/pixel coordinates) or ‘wcs’ (to assume WCS,
     usually RA/Dec, coordinates), see *note Crop modes:: for a full
     description.

‘-w FLT[,FLT[,...]]’
‘--width=FLT[,FLT[,...]]’
     Width of the cropped region about coordinate given to ‘--center’.
     If in WCS mode, value(s) given to this option will be read in the
     same units as the dataset's WCS information along this dimension
     (unless ‘--widthinpix’ is given).  This option may take either a
     single value (to be used for all dimensions: ‘--width=10’ in
     image-mode will crop a $10\times10$ pixel image) or multiple values
     (a specific value for each dimension: ‘--width=10,20’ in image-mode
     will crop a $10\times20$ pixel image).

     The ‘--width’ option also accepts fractions.  For example, if you
     want the width of your crop to be 3 by 5 arcseconds along RA and
     Dec respectively and you are in wcs-mode, you can use:
     ‘--width=3/3600,5/3600’.

     The final output will have an odd number of pixels to allow easy
     identification of the pixel which keeps your requested coordinate
     (from ‘--center’ or ‘--catalog’).  If you want an even sided crop,
     you can run Crop afterwards with ‘--section=":*-1,:*-1"’ or
     ‘--section=2:,2:’ (depending on which side you do not need), see
     *note Crop section syntax::.

     The basic reason for making an odd-sided crop is that your given
     central coordinate will ultimately fall within a discrete pixel in
     the image (defined by the FITS standard).  When the crop has an odd
     number of pixels in each dimension, that pixel can be very well
     defined as the "central" pixel of the crop, making it unambiguously
     easy to identify.  However, for an even-sided crop, it will be very
     hard to identify the central pixel (it can be on any of the four
     pixels adjacent to the central point of the image!).

‘-X’
‘--widthinpix’
     In WCS mode, interpret the value to ‘--width’ as number of pixels,
     not the WCS units like degrees.  This is useful when you want a
     fixed crop size in pixels, even though your center coordinates are
     in WCS (for example, RA and Dec).

‘-l STR’
‘-l FLT:FLT,...’
‘--polygon=STR’
‘--polygon=FLT,FLT:FLT,FLT:...’
     Polygon vertice coordinates (when value is in ‘FLT,FLT:FLT,FLT:...’
     format) or the filename of a SAO DS9 region file (when the value
     has no ‘,’ or ‘:’ characters).  Each vertice can either be in
     degrees (a single floating point number) or sexagesimal (in formats
     of '‘_h_m_’' for RA and '‘_d_m_’' for Dec, or simply '‘_:_:_’' for
     either of them).

     The vertices are used to define the polygon: in the same order
     given to this option.  When the vertices are not necessarily
     ordered in the proper order (for example, one vertice in a square
     comes after its diagonal opposite), you can add the ‘--polygonsort’
     option which will attempt to sort the vertices before cropping.
     Note that for concave polygons, sorting is not recommended because
     there is no unique solution, for more, see the description under
     ‘--polygonsort’.

     This option can be used both in the image and WCS modes, see *note
     Crop modes::.  If a SAO DS9 region file is used, the coordinate
     mode of Crop will be determined by the contents of the file and any
     value given to ‘--mode’ is ignored.  The cropped image will be the
     size of the rectangular region that completely encompasses the
     polygon.  By default all the pixels that are outside of the polygon
     will be set as blank values (see *note Blank pixels::).  However,
     if ‘--polygonout’ is called all pixels internal to the vertices
     will be set to blank.  In WCS-mode, you may provide many FITS
     images/tiles: Crop will stitch them to produce this cropped region,
     then apply the polygon.

     The syntax for the polygon vertices is similar to, and simpler
     than, that for ‘--section’.  In short, the dimensions of each
     coordinate are separated by a comma (<,>) and each vertex is
     separated by a colon (<:>).  You can define as many vertices as you
     like.  If you would like to use space characters between the
     dimensions and vertices to make them more human-readable, then you
     have to put the value to this option in double quotation marks.

     For example, let's assume you want to work on the deepest part of
     the WFC3/IR images of Hubble Space Telescope eXtreme Deep Field
     (HST-XDF). According to the web page
     (https://archive.stsci.edu/prepds/xdf/)(2) the deepest part is
     contained within the coordinates:

          [ (53.187414,-27.779152), (53.159507,-27.759633),
            (53.134517,-27.787144), (53.161906,-27.807208) ]

     They have provided mask images with only these pixels in the
     WFC3/IR images, but what if you also need to work on the same
     region in the full resolution ACS images?  Also what if you want to
     use the CANDELS data for the shallow region?  Running Crop with
     ‘--polygon’ will easily pull out this region of the image for you,
     irrespective of the resolution.  If you have set the operating mode
     to WCS mode in your nearest configuration file (see *note
     Configuration files::), there is no need to call ‘--mode=wcs’ on
     the command-line.

          $ astcrop --mode=wcs desired-filter-image(s).fits           \
             --polygon="53.187414,-27.779152 : 53.159507,-27.759633 : \
                        53.134517,-27.787144 : 53.161906,-27.807208"

     More generally, you have an image and want to define the polygon
     yourself (it is not already published like the example above).  As
     the number of vertices increases, checking the vertex coordinates
     on a FITS viewer (for example, SAO DS9) and typing them in, one by
     one, can be very tedious and prone to typo errors.  In such cases,
     you can make a polygon "region" in DS9 and using your mouse, easily
     define (and visually see) it.  Given that SAO DS9 has a graphic
     user interface (GUI), if you do not have the polygon vertices
     before-hand, it is much more easier build your polygon there and
     pass it onto Crop through the region file.

     You can take the following steps to make an SAO DS9 region file
     containing your polygon.  Open your desired FITS image with SAO DS9
     and activate its "region" mode with Edit→Region.  Then define the
     region as a polygon with Region→Shape→Polygon.  Click on the
     approximate center of the region you want and a small square will
     appear.  By clicking on the vertices of the square you can shrink
     or expand it, clicking and dragging anywhere on the edges will
     enable you to define a new vertex.  After the region has been
     nicely defined, save it as a file with Region→"Save Regions".  You
     can then select the name and address of the output file, keep the
     format as ‘REG (*.reg)’ and press the "OK" button.  In the next
     window, keep format as "ds9" and "Coordinate System" as "icrs" or
     "fk5" for RA and Dec (or "Image" for pixel coordinates).  A plain
     text file is now created (let's call it ‘ds9.reg’) which you can
     pass onto Crop with ‘--polygon=ds9.reg’.

     For the expected format of the region file, see the description of
     ‘gal_ds9_reg_read_polygon’ in *note SAO DS9 library::.  However,
     since SAO DS9 makes this file for you, you do not usually need to
     worry about its internal format unless something un-expected
     happens and you find a bug.

‘--polygonout’
     Keep all the regions outside the polygon and mask the inner ones
     with blank pixels (see *note Blank pixels::).  This is practically
     the inverse of the default mode of treating polygons.  Note that
     this option only works when you have only provided one input image.
     If multiple images are given (in WCS mode), then the full area
     covered by all the images has to be shown and the polygon excluded.
     This can lead to a very large area if large surveys like COSMOS are
     used.  So Crop will abort and notify you.  In such cases, it is
     best to crop out the larger region you want, then mask the smaller
     region with this option.

‘--polygonsort’
     Sort the given set of vertices to the ‘--polygon’ option.  For a
     concave polygon it will sort the vertices correctly, however for a
     convex polygon it there is no unique sorting, so be careful because
     the crop may not be what you expected.

     Polygons come in two classes: convex and concave (or generally,
     non-convex!), see below for a demonstration.  Convex polygons are
     those where all inner angles are less than 180 degrees.  By
     contrast, a concave polygon is one where an inner angle may be more
     than 180 degrees.

                      Concave Polygon        Convex Polygon

                       D --------C          D------------- C
                        \        |        E /              |
                         \E      |          \              |
                         /       |           \             |
                        A--------B             A ----------B

‘-s STR’
‘--section=STR’
     Section of the input image which you want to be cropped.  See *note
     Crop section syntax:: for a complete explanation on the syntax
     required for this input.

‘-C FITS/TXT’
‘--catalog=FITS/TXT’
     File name of catalog for making multiple crops from the input
     images/cubes.  The catalog can be in any of Gnuastro's recognized
     *note Recognized table formats::.  The columns containing the
     coordinates for the crop centers can be specified with the
     ‘--coordcol’ option (using column names or numbers, see *note
     Selecting table columns::).  The catalog can also contain the name
     of each crop, you can specify the column containing the name with
     the ‘--namecol’.

‘--cathdu=STR/INT’
     The HDU (extension) containing the catalog (if the file given to
     ‘--catalog’ is a FITS file).  This can either be the HDU name (if
     it has one) or number (counting from 0).  By default (if this
     option is not given), the second HDU will be used (equivalent to
     ‘--cathdu=1’.  For more on how to specify the HDU, see the
     explanation of the ‘--hdu’ option in *note Input output options::.

‘-x STR/INT’
‘--coordcol=STR/INT’
     The column in a catalog to read as a coordinate.  The value can be
     either the column number (starting from 1), or a match/search in
     the table meta-data, see *note Selecting table columns::.  This
     option must be called multiple times, depending on the number of
     dimensions in the input dataset.  If it is called more than
     necessary, the extra columns (later calls to this option on the
     command-line or configuration files) will be ignored, see *note
     Configuration file precedence::.

‘-n STR/INT’
‘--namecol=STR/INT’
     Column selection of crop file name.  The value can be either the
     column number (starting from 1), or a match/search in the table
     meta-data, see *note Selecting table columns::.  This option can be
     used both in Image and WCS modes, and not a mandatory.  When a
     column is given to this option, the final crop base file name will
     be taken from the contents of this column.  The directory will be
     determined by the ‘--output’ option (current directory if not
     given) and the value to ‘--suffix’ will be appended.  When this
     column is not given, the row number will be used instead.

   ---------- Footnotes ----------

   (1) <https://archive.stsci.edu/pub/hlsp/candels/goods-s/gs-tot/v1.0/>

   (2) <https://archive.stsci.edu/prepds/xdf/>


File: gnuastro.info,  Node: Crop output,  Next: Crop known issues,  Prev: Crop options,  Up: Invoking astcrop

6.1.4.2 Crop output
...................

The string given to ‘--output’ option will be interpreted depending on
how many crops were requested, see *note Crop modes:::

   • When a catalog is given, the value of the ‘--output’ (see *note
     Common options::) will be read as the directory to store the output
     cropped images.  Hence if it does not already exist, Crop will
     abort with an "No such file or directory" error.

     The crop file names will consist of two parts: a variable part (the
     row number of each target starting from 1) along with a fixed
     string which you can set with the ‘--suffix’ option.  Optionally,
     you may also use the ‘--namecol’ option to define a column in the
     input catalog to use as the file name instead of numbers.

   • When only one crop is desired, the value to ‘--output’ will be read
     as a file name.  If no output is specified or if it is a directory,
     the output file name will follow the automatic output names of
     Gnuastro, see *note Automatic output::: The string given to
     ‘--suffix’ will be replaced with the ‘.fits’ suffix of the input.

   When the desired crop is not within the input image(s) crop will not
produce any image for that crop.  If the ‘--quiet’ option is not given,
crop will report which one of the inputs was created and which wasn't
(due to an overlap or the center being empty, see ‘--checkcenter’.  This
information can be formally written in an optional log-file also, see
below.  At the end, Crop will return successfully to the shell if at
least one output file was created.  If no output was created at all,
then crop will return to the shell with a failure.

   By default, as suggested by the FITS standard and implemented in all
Gnuastro programs, the first/primary extension of the output files will
only contain metadata.  The cropped images/cubes will be written into
the 2nd HDU of their respective FITS file (which is actually counted as
‘1’ because HDU counting starts from ‘0’).  However, if you want the
cropped data to be written into the primary (0-th) HDU, run Crop with
the ‘--primaryimghdu’ option.

   If the output file already exists by default Crop will re-write it
(so that all existing HDUs in it will be deleted).  If you want the
cropped HDU to be appended to existing HDUs, use ‘--append’ described
below.

   The 0-th HDU of each output cropped image will contain the names of
the input image(s) it was cut from.  If a name is longer than the 70
character space that the FITS standard allows for header keyword values,
the name will be cut into several keywords from the nearest slash (</>).
The keywords have the following format: ‘ICFn_m’ (for Crop File).  Where
‘n’ is the number of the image used in this crop and ‘m’ is the part of
the name (it can be broken into multiple keywords).  Following the name
is another keyword named ‘ICFnPIX’ which shows the pixel range from that
input image in the same syntax as *note Crop section syntax::.  So this
string can be directly given to the ‘--section’ option later.

   Once done, a log file can be created in the current directory with
the ‘--log’ option.  This file will have three columns and the same
number of rows as the number of cropped images.  There are also comments
on the top of the log file explaining basic information about the run
and descriptions for the columns.  A short description of the columns is
also given below:

  1. The cropped image file name for that row.
  2. The number of input images that were used to create that image.
  3. A ‘0’ if the central few pixels (value to the ‘--checkcenter’
     option) are blank and ‘1’ if they are not.  When the crop was not
     defined by its center (see *note Crop modes::), or ‘--checkcenter’
     was given a value of 0 (see *note Invoking astcrop::), the center
     will not be checked and this column will be given a value of ‘-1’.

   If the output crop(s) have a single element (pixel in an image) and
‘--oneelemstdout’ has been called, no output file will be produced!
Instead, the single element's value is printed on the standard output.
See the description of ‘--oneelemstdout’ below for more:

‘-p STR’
‘--suffix=STR’
     The suffix (or post-fix) of the output files for when you want all
     the cropped images to have a special ending.  One case where this
     might be helpful is when besides the science images, you want the
     weight images (or exposure maps, which are also distributed with
     survey images) of the cropped regions too.  So in one run, you can
     set the input images to the science images and ‘--suffix=_s.fits’.
     In the next run you can set the weight images as input and
     ‘--suffix=_w.fits’.

‘-a STR’
‘--metaname=STR’
     Name of cropped HDU (value to the ‘EXTNAME’ keyword of FITS). If
     not given, a default ‘CROP’ will be placed there (so the ‘EXTNAME’
     keyword will always be present in the output).  If crop produces
     many outputs from a catalog, they will be given the same string as
     ‘EXTNAME’ (the file names containing the cropped HDU will be
     different).

‘-A’
‘--append’
     If the output file already exists, append the cropped image HDU to
     the end of any existing HDUs.  By default (when this option isn't
     given), if an output file already exists, any existing HDU in it
     will be deleted.  If the output file doesn't exist, this option is
     redundant.

‘--primaryimghdu’
     Write the output into the primary (0-th) HDU/extension of the
     output.  By default, like all Gnuastro's default outputs, no data
     is written in the primary extension because the FITS standard
     suggests keeping that extension free of data and only for metadata.

‘-t’
‘--oneelemstdout’
     When a crop only has a single element (a single pixel), print it to
     the standard output instead of making a file.  By default (without
     this option), a single-pixel crop will be saved to a file, just
     like a crop of any other size.

     When a single crop is requested (either through ‘--center’, or a
     catalog of one row is given), the single value alone is printed
     with nothing else.  This makes it easy to immediately write the
     value into a shell variable for example:

          value=$(astcrop img.fits --mode=wcs --center=1.234,5.678 \
                          --width=1 --widthinpix --oneelemstdout \
                          --quiet)

     If a catalog of coordinates is given (that would produce multiple
     crops; or multiple values in this scenario), the solution for a
     single value will not work!  Recall that Crop will do the crops in
     parallel, therefore each time you run it, the order of the rows
     will be different and not correspond to the order of the inputs.

     To allow identification of each value (which row of the input
     catalog it corresponds to), Crop will first print the name of the
     would-be created file name, and print the value after it (separated
     by an empty SPACE character).  In other words, the file in the
     first column will not actually be created, but the value of the
     pixel it would have contained (if this option was not called) is
     printed after it.

‘-c FLT/INT’
‘--checkcenter=FLT/INT’
     Square box width of region in the center of the image to check for
     blank values.  If any of the pixels in this central region of a
     crop (defined by its center) are blank, then it will not be stored
     in an output file.  If the value to this option is zero, no
     checking is done.  This check is only applied when the cropped
     region(s) are defined by their center (not by the vertices, see
     *note Crop modes::) and when FITS files are made (‘--oneelemstdout’
     is not called).

     The units of the value are interpreted based on the ‘--mode’ value
     (in WCS or pixel units).  The ultimate checked region size (in
     pixels) will be an odd integer around the center (converted from
     WCS, or when an even number of pixels are given to this option).
     In WCS mode, the value can be given as fractions, for example, if
     the WCS units are in degrees, ‘0.1/3600’ will correspond to a check
     size of 0.1 arcseconds.

     Because survey regions do not often have a clean square or
     rectangle shape, some of the pixels on the sides of the survey FITS
     image do not commonly have any data and are blank (see *note Blank
     pixels::).  So when the catalog was not generated from the input
     image, it often happens that the image does not have data over some
     of the points.

     When the given center of a crop falls in such regions or outside
     the dataset, and this option has a non-zero value, no crop will be
     created.  Therefore with this option, you can specify a width of a
     small box (3 pixels is often good enough) around the central pixel
     of the cropped image.  You can check which crops were created and
     which were not from the command-line (if ‘--quiet’ was not called,
     see *note Operating mode options::), or in Crop's log file (see
     *note Crop output::).

‘-b’
‘--noblank’
     Pixels outside of the input image that are in the crop box will not
     be used.  By default they are filled with blank values (depending
     on type), see *note Blank pixels::.  This option only applies only
     in Image mode, see *note Crop modes::.

‘-z’
‘--zeroisnotblank’
     In float or double images, it is common to give the value of zero
     to blank pixels.  If the input image type is one of these two
     types, such pixels will also be considered as blank.  You can
     disable this behavior with this option, see *note Blank pixels::.

‘--log’
     Generate a log file that contains the status of each of each crop
     in three columns: the name of the crop, if its center overlapped
     with the input(s) and if the central pixels were blank or not.  The
     name of the log file will be based on the output file name:
        • If no output file name is given, it will be
          ‘astcrop-log.fits’.
        • If there is only a single crop and the output is a FITS image,
          the log file will have the same base name as the output crop,
          but with a ‘-log.fits’ suffix.  In case ‘--oneelemstdout’ is
          called with a single-pixel width then no output file will be
          created, so the log-file's name will be the string given to
          ‘--output’.
        • If the output is a directory name (when a catalog is given),
          the log file name will be the directory name along with a
          ‘-log.fits’ suffix.  In case ‘--oneelemstdout’ is called with
          a single-pixel width then no output file will be created, so
          the log-file's name will be the string given to ‘--output’.


File: gnuastro.info,  Node: Crop known issues,  Prev: Crop output,  Up: Invoking astcrop

6.1.4.3 Crop known issues
.........................

When running Crop, you may encounter strange errors and bugs.  In these
cases, please report a bug and we will try to fix it as soon as
possible, see *note Report a bug::.  However, some things are beyond our
control, or may take too long to fix directly.  In this section we list
such known issues that may occur in known cases and suggest the hack (or
work-around) to fix the problem:

Crash with ‘Killed’ when cropping catalog from ‘.fits.gz’
     This happens because CFISTIO (that reads and writes FITS files)
     will internally decompress the file in a temporary place (possibly
     in the RAM), then start reading from it.  On the other hand, by
     default when given a catalog (with many crops) and not specifying
     ‘--numthreads’, Crop will use the maximum number of threads
     available on your system to do each crop faster.  On an normal (not
     compressed) file, parallel access will not cause a problem,
     however, when attempting parallel access with the maximum number of
     threads on a compressed file, CFITSIO crashes with ‘Killed’.
     Therefore the following solutions can be used to fix this crash:

        • Decrease the number of threads (at the minimum, set
          ‘--numthreads=1’).  Since this solution does not attempt to
          change any of your previous Crop command components or does
          not change your local file structure, it is the preferred way.

        • Decompress the file (with the command below) and feed the
          ‘.fits’ file into Crop without changing the number of threads.

               $ gunzip -k image.fits.gz


File: gnuastro.info,  Node: Arithmetic,  Next: Convolve,  Prev: Crop,  Up: Data manipulation

6.2 Arithmetic
==============

It is commonly necessary to do operations on some or all of the elements
of a dataset independently (pixels in an image).  For example, in the
reduction of raw data it is necessary to subtract the Sky value (*note
Sky value::) from each image image.  Later (once the images as warped
into a single grid using Warp for example, see *note Warp::), the images
are co-added (the output pixel grid is the average of the pixels of the
individual input images).  Arithmetic is Gnuastro's program for such
operations on your datasets directly from the command-line.  It
currently uses the reverse polish or post-fix notation, see *note
Reverse polish notation:: and will work on the native data types of the
input images/data to reduce CPU and RAM resources, see *note Numeric
data types::.  For more information on how to run Arithmetic, please see
*note Invoking astarithmetic::.

* Menu:

* Reverse polish notation::     The current notation style for Arithmetic.
* Integer benefits and pitfalls::  Integers have benefits, but require care.
* Noise basics::                Introduction various noise models.
* Arithmetic operators::        List of operators known to Arithmetic.
* Invoking astarithmetic::      How to run Arithmetic: options and output.


File: gnuastro.info,  Node: Reverse polish notation,  Next: Integer benefits and pitfalls,  Prev: Arithmetic,  Up: Arithmetic

6.2.1 Reverse polish notation
-----------------------------

The most common notation for arithmetic operations is the infix notation
(https://en.wikipedia.org/wiki/Infix_notation) where the operator goes
between the two operands, for example, $4+5$.  The infix notation is the
preferred way in most programming languages which come with scripting
features for large programs.  This is because the infix notation
requires a way to define precedence when more than one operator is
involved.

   For example, consider the statement ‘5 + 6 / 2’.  Should 6 first be
divided by 2, then added by 5?  Or should 5 first be added with 6, then
divided by 2?  Therefore we need parenthesis to show precedence:
‘5+(6/2)’ or ‘(5+6)/2’.  Furthermore, if you need to leave a value for
later processing, you will need to define a variable for it; for
example, ‘a=(5+6)/2’.

   Gnuastro provides libraries where you can also use infix notation in
C or C++ programs.  However, Gnuastro's programs are primarily designed
to be run on the command-line and the level of complexity that infix
notation requires can be annoying/confusing to write on the command-line
(where they can get confused with the shell's parenthesis or variable
definitions).  Therefore Gnuastro's Arithmetic and Table (when doing
column arithmetic) programs use the post-fix notation, also known as
reverse polish notation
(https://en.wikipedia.org/wiki/Reverse_Polish_notation).  For example,
instead of writing ‘5+6’, we write ‘5 6 +’.

   The Wikipedia article on the reverse polish notation provides some
excellent explanation on this notation but here we will give a short
summary here for self-sufficiency.  In short, in the reverse polish
notation, the operator is placed after the operands.  As we will see
below this removes the need to define parenthesis and lets you use
previous values without needing to define a variable.  In the future(1)
we do plan to also optionally allow infix notation when arithmetic
operations on datasets are desired, but due to time constraints on the
developers we cannot do it immediately.

   To easily understand how the reverse polish notation works, you can
think of each operand (‘5’ and ‘6’ in the example above) as a node in a
"last-in-first-out" stack.  One such stack in daily life is a stack of
dishes in the kitchen: you put a clean dish, on the top of a stack of
dishes when it is ready for later usage.  Later, when you need a dish,
you pick the top one (hence the "last" dish placed "in" the stack is the
"first" dish that comes "out" when necessary).

   Each operator will need a certain number of operands (in the example
above, the ‘+’ operator needs two operands: ‘5’ and ‘6’).  In the
kitchen metaphor, an operator can be an oven.  Every time an operator is
confronted, the operator takes (or "pops") the number of operands it
needs from the top of the stack (so they do not exist in the stack any
more), does its operation, and places (or "pushes") the result back on
top of the stack.  So if you want the average of 5 and 6, you would
write: ‘5 6 + 2 /’.  The operations that are done are:

  1. ‘5’ is an operand, so Arithmetic pushes it to the top of the stack
     (which is initially empty).  In the kitchen metaphor, you can
     visualize this as taking a new dish from the cabinet, putting the
     number 5 inside of the dish, and putting the dish on top of the
     (empty) cooking table in front of you.  You now have a stack of one
     dish on the table in front of you.
  2. ‘6’ is also an operand, so it is pushed to the top of the stack.
     Like before, you can visualize this as taking a new dish from the
     cabinet, putting the number 6 in it and placing it on top of the
     previous dish.  You now have a stack of two dishes on the table in
     front of you.
  3. ‘+’ is a _binary_ operator, so it will pop the top two elements of
     the stack out of it, and perform addition on them (the order is
     $5+6$ in the example above).  The result is ‘11’ which is pushed to
     the top of the stack.

     To visualize this, you can think of the ‘+’ operator as an oven
     with a place for two dishes.  You pick up the top-most dish (that
     has the number 6 in it) and put it in the oven.  The top dish is
     now the one that has the number 5.  You also pick it up and put it
     in the oven, and close the oven door.  When the oven has finished
     its cooking, it produces a single output (in one dish, with the
     number 11 inside of it).  You take that output dish and put it back
     on the table.  You now have a stack of one dish on the table in
     front of you.
  4. ‘2’ is an operand so push it onto the top of the stack.  In the
     kitchen metaphor, you again go to the cabinet, pick up a dish and
     put the number 2 inside of it and put the dish over the previous
     dish (that has the number 11).  You now have a stack of two dishes
     on the table in front of you.
  5. ‘/’ (division) is a binary operator, so pull out the top two
     elements of the stack (top-most is ‘2’, then ‘11’) and divide the
     second one by the first.  In the kitchen metaphor, the ‘/’ operator
     can be visualized as a microwave that takes two dishes.  But unlike
     the oven (‘+’ operator) before, the order of inputs matters (they
     are on top of each other: with the top dish holder being the
     numerator and the bottom one being the denominator).  Again, you
     look at your stack of dishes on the table.

     You pick up the top one (with value 2 inside of it) and put it in
     the microwave's bottom (denominator) dish holder.  Then you go back
     to your stack of dishes on the table and pick up the top dish (with
     value 11 inside of it) and put that in the top (nominator) dish
     holder.  The microwave will do its work and when it is finished,
     returns a new dish with the single value 5.5 inside of it.  You
     pick up the dish from the microwave and place it back on the table.

  6. There are no more operands or operators, so simply return the
     remaining operand in the output.  In the kitchen metaphor, you see
     that your recipe has no more steps, so you just pick up the
     remaining dish and take it to the dining room to enjoy a good
     dinner.

   In the Arithmetic program, the operands can be FITS images of any
dimensionality, or numbers (see *note Invoking astarithmetic::).  In
Table's column arithmetic, they can be any column in the table (a series
of numbers in an array) or a single number (see *note Column
arithmetic::).

   With this notation, very complicated procedures can be created
without the need for parenthesis or worrying about precedence.  Even
functions which take an arbitrary number of arguments can be defined in
this notation.  This is a very powerful notation and is used in
languages like Postscript (2) which produces PDF files when compiled.

   ---------- Footnotes ----------

   (1) <https://savannah.gnu.org/task/index.php?13867>

   (2) See the EPS and PDF part of *note Recognized file formats:: for a
little more on the Postscript language.


File: gnuastro.info,  Node: Integer benefits and pitfalls,  Next: Noise basics,  Prev: Reverse polish notation,  Up: Arithmetic

6.2.2 Integer benefits and pitfalls
-----------------------------------

Integers are the simplest numerical data types (*note Numeric data
types::).  Because of this, their storage space is much less, and their
processing is much faster than floating point types.  You can confirm
this on your computer with the series of commands below.  You will make
four 5000 by 5000 pixel images filled with random values.  Two of them
will be saved as signed 8-bit integers, and two with 64-bit floating
point types.  The last command prints the size of the created images.

     $ astarithmetic 5000 5000 2 makenew 5 mknoise-sigma int8    -oint-1.fits
     $ astarithmetic 5000 5000 2 makenew 5 mknoise-sigma int8    -oint-2.fits
     $ astarithmetic 5000 5000 2 makenew 5 mknoise-sigma float64 -oflt-1.fits
     $ astarithmetic 5000 5000 2 makenew 5 mknoise-sigma float64 -oflt-2.fits
     $ ls -lh int-*.fits flt-*.fits

   The 8-bit integer images are only 24MB, while the 64-bit floating
point images are 191 MB! Besides helping in storage (on your disk, or in
RAM, while the program is running), the small size of these files also
helps in faster reading of the inputs.  Furthermore, CPUs can process
integer operations much faster than floating points.  In the integers,
the ones with a smaller width (number of bits) can be processed much
faster.  You can see this with the two commands below where you will add
the integer images with each other and the floats with each other:

     $ astarithmetic flt-1.fits flt-2.fits + -oflt-sum.fits -g1
     $ astarithmetic int-1.fits int-2.fits + -oint-sum.fits -g1

   Have a look at the running time of the two commands above (that is
printed on their last line).  On the system that this paragraph was
written on, the floating point and integer image sums were respectively
done in 0.481 and 0.089 seconds (the integer operation was almost 5
times faster!).

*If your data does not have decimal points, use integer types:* integer
types are much faster and can take much less space in your storage or
RAM (while the program is running).

*Select the smallest width that can host the range/precision of values:*
for example, if the largest possible value in your dataset is 1000 and
all numbers are integers, store it as a 16-bit integer.  Also, if you
know the values can never become negative, store it as an unsigned
16-bit integer.  For floating point types, if you know you will not need
a precision of more than 6 significant digits, use the 32-bit floating
point type.  For more on the range (for integers) and precision (for
floats), see *note Numeric data types::.

   There is a price to be paid for this improved efficiency in integers:
your wisdom!  If you have not selected your types wisely, strange
situations may happen.  For example, try the command below:

     $ astarithmetic 125 10 +

You expect the output to be $135$, but it will be $-121$!  The reason is
that when Arithmetic (or column-arithmetic in Table) confronts a number
on the command-line, it use the principles above to select the most
efficient type for each number.  Both $125$ and $10$ can safely fit
within a signed, 8-bit integer type, so arithmetic will store both as an
8-bit integer.  However, the sum ($135$) is larger than the maximum
possible value of an 8-bit signed integer ($127$).  Therefore an integer
overflow will occur, and the bits will be over-written.  As a result,
the value will be $135-128=7$ more than the minimum value of this type
($-128$), which is $-128+7=-121$.

   When you know situations like this may occur, you can simply use
*note Numerical type conversion operators::, to set just one of the
inputs to a wider data type (the smallest, wider type to avoid wasting
resources).  In the example above, this would be ‘uint16’:

     $ astarithmetic 125 uint16 10 +

   The reason this worked is that $125$ is now converted into an
unsigned 16-bit integer before the ‘+’ operator.  Since this is larger
than an 8-bit integer, the C programming language's automatic type
conversion will treat both as the wider type and store the result of the
binary operation (‘+’) in that type.

   For such a basic operation like the command above, a faster hack
would be any of the two commands below (which are equivalent).  This is
because ‘125.0’ or ‘125.’ are interpreted as floating-point types and
they do not suffer from such issues (converting only on one input is
enough):

     $ astarithmetic 125.  10 +
     $ astarithmetic 125.0 10 +

   For this particular command, the fix above will be as fast as the
‘uint16’ solution.  This is because there are only two numbers, and the
overhead of Arithmetic (reading configuration files, etc.)  dominates
the running time.  However, for large datasets, the ‘uint16’ solution
will be faster (as you saw above), Arithmetic will consume less RAM
while running, and the output will consume less storage in your system
(all major benefits)!

   It is possible to do internal checks in Gnuastro and catch integer
overflows and correct them internally.  However, we have not opted for
this solution because all those checks will consume significant
resources and slow down the program (especially with large datasets
where RAM, storage and running time become important).  To be optimal,
we therefore trust that you (the wise Gnuastro user!)  make the
appropriate type conversion in your commands where necessary (recall
that the operators are available in *note Numerical type conversion
operators::).


File: gnuastro.info,  Node: Noise basics,  Next: Arithmetic operators,  Prev: Integer benefits and pitfalls,  Up: Arithmetic

6.2.3 Noise basics
------------------

Deep astronomical images, like those used in extragalactic studies,
seriously suffer from noise in the data.  Generally speaking, the
sources of noise in an astronomical image are photon counting noise and
Instrumental noise which are discussed in *note Photon counting noise::
and *note Instrumental noise::.  This review finishes with *note
Generating random numbers:: which is a short introduction on how random
numbers are generated.  We will see that while software random number
generators are not perfect, they allow us to obtain a reproducible
series of random numbers through setting the random number generator
function and seed value.  Therefore in this section, we will also
discuss how you can set these two parameters in Gnuastro's programs
(including the arithmetic operators in *note Random number
generators::).

* Menu:

* Photon counting noise::       Poisson noise
* Instrumental noise::          Readout, dark current and other sources.
* Final noised pixel value::    How the final noised value is calculated.
* Generating random numbers::   How random numbers are generated.


File: gnuastro.info,  Node: Photon counting noise,  Next: Instrumental noise,  Prev: Noise basics,  Up: Noise basics

6.2.3.1 Photon counting noise
.............................

With the very accurate electronics used in today's detectors, photon
counting noise(1) is the most significant source of uncertainty in most
datasets.  To understand this noise (error in counting) and its effect
on the images of astronomical targets, let's start by reviewing how a
distribution produced by counting can be modeled as a parametric
function.

   Counting is an inherently discrete operation, which can only produce
positive integer outputs (including zero).  For example, we cannot count
$3.2$ or $-2$ of anything.  We only count $0$, $1$, $2$, $3$ and so on.
The distribution of values, as a result of counting efforts is formally
known as the Poisson distribution
(https://en.wikipedia.org/wiki/Poisson_distribution).  It is associated
to Siméon Denis Poisson, because he discussed it while working on the
number of wrongful convictions in court cases in his 1837 book(2).

   Let's take $\lambda$ to represent the expected mean count of
something.  Furthermore, let's take $k$ to represent the output of a
counting attempt (hence $k$ is a positive integer).  The probability
density function of getting $k$ counts (in each attempt, given the
expected/mean count of $\lambda$) can be written as:

$$f(k)={\lambda^k \over k!} e^{-\lambda},\quad k\in {0, 1, 2, 3, \dots }$$

   Because the Poisson distribution is only applicable to positive
integer values (note the factorial operator, which only applies to
non-negative integers), naturally it is very skewed when $\lambda$ is
near zero.  One qualitative way to understand this behavior is that for
smaller values near zero, there simply are not enough integers smaller
than the mean, than integers that are larger.  Therefore to accommodate
all possibilities/counts, it has to be strongly skewed to the positive
when the mean is small.  For more on Skewness, see *note Skewness caused
by signal and its measurement::.

   As $\lambda$ becomes larger, the distribution becomes more and more
symmetric, and the variance of that distribution is equal to its mean.
In other words, the standard deviation is the square root of the mean.
It can also be proved that when the mean is large, say $\lambda>1000$,
the Poisson distribution approaches the Normal (Gaussian) distribution
(https://en.wikipedia.org/wiki/Normal_distribution) with mean
$\mu=\lambda$ and standard deviation $\sigma=\sqrt{\lambda}$.  In other
words, a Poisson distribution (with a sufficiently large $\lambda$) is
simply a Gaussian that has one free parameter ($\mu=\lambda$ and
$\sigma=\sqrt{\lambda}$), instead of the two parameters that the
Gaussian distribution originally has (independent $\mu$ and $\sigma$).

   In real situations, the photons/flux from our targets are combined
with photons from a certain background (observationally, the _Sky_
value).  The Sky value is defined to be the average flux of a region in
the dataset with no targets.  Its physical origin can be the brightness
of the atmosphere (for ground-based instruments), possible stray light
within the imaging instrument, the average flux of undetected targets,
etc.  The Sky value is thus an ideal definition, because in real
datasets, what lies deep in the noise (far lower than the detection
limit) is never known(3).  To account for all of these, the sky value is
defined to be the average count/value of the undetected regions in the
image.  In a mock image/dataset, we have the luxury of setting the
background (Sky) value.

   In summary, the value in each element of the dataset (pixel in an
image) is the sum of contributions from various galaxies and stars
(after convolution by the PSF, see *note PSF::).  Let's name the
convolved sum of possibly overlapping objects in each pixel as $I_{nn}$.
$nn$ represents 'no noise'.  For now, let's assume the background ($B$)
is constant and sufficiently high for the Poisson distribution to be
approximated by a Gaussian.  Then the flux of that pixel, after adding
noise, is _a random value_ taken from a Gaussian distribution with the
following mean ($\mu$) and standard deviation ($\sigma$):

            $$\mu=B+I_{nn}, \quad \sigma=\sqrt{B+I_{nn}}$$

   In astronomical instruments, $B$ is enhanced by adding a "bias" level
to each pixel before the shutter is even opened (for the exposure to
start).  As the exposure is ongoing and photo-electrons are accumulating
from the astronomical objects, a "dark" current (due to thermal
radiation of the instrument) also builds up in the pixels.  The "dark"
current will accumulate even when the shutter is closed, but the CCD
electronics are working (hence the name "dark").  This added dark level
further enhances the mean value in a real observation compared to the
raw background value (from the atmosphere for example).

   Since this type of noise is inherent in the objects we study, it is
usually measured on the same scale as the astronomical objects, namely
the magnitude system, see *note Brightness flux magnitude::.  It is then
internally converted to the flux scale for further processing.

   The equations above clearly show the importance of the background
value and its effect on the final signal to noise ratio in each pixel of
a science image.  It is therefore, one of the most important factors in
understanding the noise (and properly simulating observations where
necessary).  An inappropriately bright background value can hide the
signal of the mock profile hide behind the noise.  In other words, a
brighter background has larger standard deviation and vice versa.  As a
result, the only necessary parameter to define photon-counting noise
over a mock image of simulated profiles is the background.  For a
complete example, see *note Sufi simulates a detection::.

   To better understand the correlation between the mean (or background)
value and the noise standard deviation, let's use an analogy.  Consider
the profile of your galaxy to be analogous to the profile of a ship that
is sailing in the sea.  The height of the ship would therefore be
analogous to the maximum flux difference between your galaxy's minimum
and maximum values.  Furthermore, let's take the depth of the sea to
represent the background value: a deeper sea, corresponds to a brighter
background.  In this analogy, the "noise" would be the height of the
waves that surround the ship: in deeper waters, the waves would also be
taller (the square root of the mean depth at the ship's position).

   If the ship is in deep waters, the height of waves are greater than
when the ship is near to the beach (at lower depths).  Therefore, when
the ship is in the middle of the sea, there are high waves that are
capable of hiding a significant part of the ship from our perspective.
This corresponds to a brighter background value in astronomical images:
the resulting noise from that brighter background can completely wash
out the signal from a fainter galaxy, star or solar system object.

   ---------- Footnotes ----------

   (1) In practice, we are actually counting the electrons that are
produced by each photon, not the actual photons.

   (2) [From Wikipedia] Poisson's result was also derived in a previous
study by Abraham de Moivre in 1711.  Therefore some people suggest it
should rightly be called the de Moivre distribution.

   (3) In a real image, a relatively large number of very faint objects
can be fully buried in the noise and never detected.  These undetected
objects will bias the background measurement to slightly larger values.
Our best approximation is thus to simply assume they are uniform, and
consider their average effect.  See Figure 1 (a.1 and a.2) and Section
2.2 in Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664).


File: gnuastro.info,  Node: Instrumental noise,  Next: Final noised pixel value,  Prev: Photon counting noise,  Up: Noise basics

6.2.3.2 Instrumental noise
..........................

While taking images with a camera, a bias current is fed to the pixels,
the variation of the value of this bias current over the pixels, also
adds to the final image noise.  Another source of noise is the readout
noise that is produced by the electronics in the detector.
Specifically, the parts that attempt to digitize the voltage produced by
the photo-electrons in the analog to digital converter.  With the
current generation of instruments, this source of noise is not as
significant as the noise due to the background Sky discussed in *note
Photon counting noise::.

   Let $C$ represent the combined standard deviation of all these
instrumental sources of noise.  When only this source of noise is
present, the noised pixel value would be a random value chosen from a
Gaussian distribution with

            $$\mu=I_{nn}, \quad \sigma=\sqrt{C^2+I_{nn}}$$

   This type of noise is independent of the signal in the dataset, it is
only determined by the instrument.  So the flux scale (and not magnitude
scale) is most commonly used for this type of noise.  In practice, this
value is usually reported in analog-to-digital units or ADUs, not flux
or electron counts.  The gain value of the device can be used to convert
between these two, see *note Brightness flux magnitude::.


File: gnuastro.info,  Node: Final noised pixel value,  Next: Generating random numbers,  Prev: Instrumental noise,  Up: Noise basics

6.2.3.3 Final noised pixel value
................................

Based on the discussions in *note Photon counting noise:: and *note
Instrumental noise::, depending on the values you specify for $B$ and
$C$ from the above, the final noised value for each pixel is a random
value chosen from a Gaussian distribution with

          $$\mu=B+I_{nn}, \quad \sigma=\sqrt{C^2+B+I_{nn}}$$


File: gnuastro.info,  Node: Generating random numbers,  Prev: Final noised pixel value,  Up: Noise basics

6.2.3.4 Generating random numbers
.................................

As discussed above, to generate noise we need to make random samples of
a particular distribution.  So it is important to understand some
general concepts regarding the generation of random numbers.  For a very
complete and nice introduction we strongly advise reading Donald Knuth's
"The art of computer programming", volume 2, chapter 3(1).  Quoting from
the GNU Scientific Library manual, "If you do not own it, you should
stop reading right now, run to the nearest bookstore, and buy it"(2)!

   Using only software, we can only produce what is called a
psuedo-random sequence of numbers.  A true random number generator is a
hardware (let's assume we have made sure it has no systematic biases),
for example, throwing dice or flipping coins (which have remained from
the ancient times).  More modern hardware methods use atmospheric noise,
thermal noise or other types of external electromagnetic or quantum
phenomena.  All pseudo-random number generators (software) require a
seed to be the basis of the generation.  The advantage of having a seed
is that if you specify the same seed for multiple runs, you will get an
identical sequence of random numbers which allows you to reproduce the
same final noised image.

   The programs in GNU Astronomy Utilities (for example, MakeNoise or
MakeProfiles) use the GNU Scientific Library (GSL) to generate random
numbers.  GSL allows the user to set the random number generator through
environment variables, see *note Installation directory:: for an
introduction to environment variables.  In the chapter titled "Random
Number Generation" they have fully explained the various random number
generators that are available (there are a lot of them!).  Through the
two environment variables ‘GSL_RNG_TYPE’ and ‘GSL_RNG_SEED’ you can
specify the generator and its seed respectively.

   If you do not specify a value for ‘GSL_RNG_TYPE’, GSL will use its
default random number generator type.  The default type is sufficient
for most general applications.  If no value is given for the
‘GSL_RNG_SEED’ environment variable and you have asked Gnuastro to read
the seed from the environment (through the ‘--envseed’ option), then GSL
will use the default value of each generator to give identical outputs.
If you do not explicitly tell Gnuastro programs to read the seed value
from the environment variable, then they will use the system time
(accurate to within a microsecond) to generate (apparently random)
seeds.  In this manner, every time you run the program, you will get a
different random number distribution.

   There are two ways you can specify values for these environment
variables.  You can call them on the same command-line for example:

     $ GSL_RNG_TYPE="taus" GSL_RNG_SEED=345 astarithmetic input.fits \
                                                          mknoise-sigma \
                                                          --envseed

In this manner the values will only be used for this particular
execution of Arithmetic.  However, it makes your code hard to read!
Alternatively, you can define them for the full period of your terminal
session or script, using the shell's ‘export’ command with the two
separate commands below (for a script remove the ‘$’ signs):

     $ export GSL_RNG_TYPE="taus"
     $ export GSL_RNG_SEED=345

The subsequent programs which use GSL's random number generators will
hence forth use these values in this session of the terminal you are
running or while executing this script.  In case you want to set fixed
values for these parameters every time you use the GSL random number
generator, you can add these two lines to your ‘.bashrc’ startup
script(3), see *note Installation directory::.

   *IMPORTANT NOTE:* If the two environment variables ‘GSL_RNG_TYPE’ and
‘GSL_RNG_SEED’ are defined, GSL will report them by default, even if you
do not use the ‘--envseed’ option.  For example, see this call to
MakeProfiles:

     $ export GSL_RNG_TYPE=taus
     $ export GSL_RNG_SEED=345
     $ astmkprof -s1 --kernel=gaussian,2,5
     GSL_RNG_TYPE=taus
     GSL_RNG_SEED=345
     MakeProfiles V.VV started on DDD MMM DDD HH:MM:SS YYYY
       - Building one gaussian kernel
       - Random number generator (RNG) type: taus
       - Basic RNG seed: 1618960836
       ---- ./kernel.fits created.
       -- Output: ./kernel.fits
     MakeProfiles finished in 0.068945 seconds

The first two output lines (showing the names and values of the GSL
environment variables) are printed by GSL before MakeProfiles actually
starts generating random numbers.  Gnuastro's programs will report the
actual values they use independently (after the name of the program),
you should check them for the final values used, not GSL's printed
values.  In the example above, did you notice how the random number
generator seed above is different between GSL and MakeProfiles?
However, if ‘--envseed’ was given, both printed seeds would be the same.

   ---------- Footnotes ----------

   (1) Knuth, Donald.  1998.  The art of computer programming.
Addison-Wesley.  ISBN 0-201-89684-2

   (2) For students, running to the library might be more affordable!

   (3) Do not forget that if you are going to give your scripts (that
use the GSL random number generator) to others you have to make sure you
also tell them to set these environment variable separately.  So for
scripts, it is best to keep all such variable definitions within the
script, even if they are within your ‘.bashrc’.


File: gnuastro.info,  Node: Arithmetic operators,  Next: Invoking astarithmetic,  Prev: Noise basics,  Up: Arithmetic

6.2.4 Arithmetic operators
--------------------------

In this section, list of recognized operators in Arithmetic (and the
Table program's *note Column arithmetic::) and discussed in detail with
examples.  As mentioned before, to be able to easily do complex
operations on the command-line, the Reverse Polish Notation is used
(where you write '$4\quad5\quad+$' instead of '$4 + 5$'), if you are not
already familiar with it, before continuing, please see *note Reverse
polish notation::.

   The operands to all operators can be a data array (for example, a
FITS image or data cube) or a number, the output will be an array or
number according to the inputs.  For example, a number multiplied by an
array will produce an array.  The numerical data type of the output of
each operator is described within it.  Here are some generic tips and
tricks (relevant to all operators):

Multiple operators in one command
     When you need to use arithmetic commands in several consecutive
     operations, you can use one command instead of multiple commands
     and perform all calculations in the same command.  For example,
     assume you want to apply a threshold of 10 on your image, and label
     the connected groups of pixel above this threshold.  You need two
     operators for this: ‘gt’ (for "greater than", see *note Conditional
     operators::) and ‘connected-components’ (see *note Mathematical
     morphology operators::).  The bad (non-optimized and slow) way of
     doing this is to call Arithmetic two times:
          $ astarithmetic image.fits 10 gt --output=thresh.fits
          $ astarithmetic thresh.fits 2 connected-components \
                          --output=labeled.fits
          $ rm thresh.fits

     The good (optimal) way is to call them after each other (remember
     *note Reverse polish notation::):
          $ astarithmetic image.fits 10 gt 2 connected-components \
                          --output=labeled.fits

     You can similarly add any number of operations that must be done
     sequentially in a single command and benefit from the speed and
     lack of intermediate files.  When your commands become long, you
     can use the ‘set-AAA’ operator to make it more readable, see *note
     Operand storage in memory or a file::.

Blank pixels in Arithmetic
     Blank pixels in the image (see *note Blank pixels::) will be stored
     based on the data type.  When the input is floating point type,
     blank values are NaN. One aspect of NaN values is that by
     definition they will fail on _any_ comparison.  Also, any operator
     that includes a NaN as a an operand will produce a NaN
     (irrespective of its other operands).  Hence both equal and
     not-equal operators will fail when both their operands are NaN!
     Therefore, the only way to guarantee selection of blank pixels is
     through the ‘isblank’ operator explained above.

     One way you can exploit this property of the NaN value to your
     advantage is when you want a fully zero-valued image (even over the
     blank pixels) based on an already existing image (with same size
     and world coordinate system settings).  The following command will
     produce this for you:

          $ astarithmetic input.fits nan eq --output=all-zeros.fits

     Note that on the command-line you can write NaN in any case (for
     example, ‘NaN’, or ‘NAN’ are also acceptable).  Reading NaN as a
     floating point number in Gnuastro is not case-sensitive.

* Menu:

* Basic mathematical operators::  For example, +, -, /, log, and pow.
* Trigonometric and hyperbolic operators::  sin, cos, atan, asinh, etc.
* Constants::                   Physical and Mathematical constants.
* Coordinate conversion operators::  For example equatorial J2000 to Galactic.
* Unit conversion operators::   Various unit conversions necessary.
* Statistical operators::       Statistics of a single dataset (for example, mean).
* Coadding operators::          Coadding or combining multiple datasets into one.
* Filtering operators::         Smoothing a dataset through mixing pixel with neighbors.
* Pooling operators::           Reducing size through statistics of pixels in window.
* Interpolation operators::     Giving blank pixels a value.
* Dimensionality changing operators::  Collapse or expand a dataset.
* Conditional operators::       Select certain pixels within the dataset.
* Mathematical morphology operators::  Work on binary images, for example, erode.
* Bitwise operators::           Work on bits within one pixel.
* Numerical type conversion operators::  Convert the numeric datatype of a dataset.
* Random number generators::    Random numbers can be used to add noise for example.
* Coordinate and border operators::  Return edges of 2D boxes.
* Loading external columns::    Read a column from a table into the stack.
* Size and position operators::  Extracting image size and pixel positions.
* New operands::                How to construct an empty dataset from scratch.
* Operand storage in memory or a file::  Tools for complex operations in one command.


File: gnuastro.info,  Node: Basic mathematical operators,  Next: Trigonometric and hyperbolic operators,  Prev: Arithmetic operators,  Up: Arithmetic operators

6.2.4.1 Basic mathematical operators
....................................

These are some of the most common operations you will be doing on your
data and include, so no further explanation is necessary.  If you are
new to Gnuastro, just read the description of each carefully.

‘+’
     Addition, so "‘4 5 +’" is equivalent to $4+5$.  For example, in the
     command below, the value 20000 is added to each pixel's value in
     ‘image.fits’:
          $ astarithmetic 20000 image.fits +
     You can also use this operator to sum the values of one pixel in
     two images (which have to be the same size).  For example, in the
     commands below (which are identical, see paragraph after the
     commands), each pixel of ‘sum.fits’ is the sum of the same pixel's
     values in ‘a.fits’ and ‘b.fits’.
          $ astarithmetic a.fits b.fits + -h1 -h1 --output=sum.fits
          $ astarithmetic a.fits b.fits + -g1     --output=sum.fits
     The HDU/extension has to be specified for each image with ‘-h’.
     However, if the HDUs are the same in all inputs, you can use ‘-g’
     to only specify the HDU once

     If you need to add more than one dataset, one way is to use this
     operator multiple times, for example, see the two commands below
     that are identical in the Reverse Polish Notation (*note Reverse
     polish notation::):
          $ astarithmetic a.fits b.fits + c.fits + -osum.fits
          $ astarithmetic a.fits b.fits c.fits + + -osum.fits

     However, this can get annoying/buggy if you have more than three or
     four images, in that case, a better way to sum data is to use the
     ‘sum’ operator (which also ignores blank pixels), that is discussed
     in *note Coadding operators::.

     *NaN values:* if a single argument of ‘+’ has a NaN value, the
     output will also be NaN. To ignore NaN values, use the ‘sum’
     operator of *note Coadding operators::.  You can see the difference
     with the two commands below:

          $ astarithmetic --quiet 1.0 2.0 3.0 nan + + +
          nan
          $ astarithmetic --quiet 1.0 2.0 3.0 nan 4 sum
          6.000000e+00

     The same goes for all the *note Coadding operators:: so if your
     data may include NaN pixels, be sure to use the coadding operators.

‘-’
     Subtraction, so "‘4 5 -’" is equivalent to $4-5$.  Usage of this
     operator is similar to ‘+’ operator, for example:
          $ astarithmetic 20000 image.fits -
          $ astarithmetic a.fits b.fits - -g1 --output=sub.fits

‘x’
     Multiplication, so "‘4 5 x’" is equivalent to $4\times5$.  For
     example, in the command below, the value of each output pixel is 5
     times its value in ‘image.fits’:
          $ astarithmetic image.fits 5 x
     And you can multiply the value of each pixel in two images, like
     this:
          $ astarithmetic a.fits a.fits x -g1 –output=multip.fits

‘/’
     Division, so "‘4 5 /’" is equivalent to $4/5$.  Like the
     multiplication, for example
          $ astarithmetic image.fits 5 -h1 /
          $ astarithmetic a.fits b.fits / -g1 –output=div.fits

‘%’
     Modulo (remainder), so "‘3 2 %’" will return $1$.  Note that the
     modulo operator only works on integer types (see *note Numeric data
     types::).  This operator is therefore not defined for most
     processed astronomical astronomical images that have floating-point
     value.  However it is useful in labeled images, for example, *note
     Segment output::).  In such cases, each pixel is the integer label
     of the object it is associated with hence with the example command
     below, we can change the labels to only be between 1 and 4 and
     decrease all objects on the image to 4/5th (all objects with a
     label that is a multiple of 5 will be set to 0).
          $ astarithmetic label.fits 5 1 %

‘abs’
     Absolute value of first operand, so "‘4 abs’" is equivalent to
     $|4|$.  For example, the output of the command bellow will not have
     any negative pixels (all negative pixels will be multiplied by $-1$
     to become positive)
          $ astarithmetic image.fits abs

‘pow’
     First operand to the power of the second, so "‘4.3 5 pow’" is
     equivalent to $4.3^{5}$.  For example, with the command below all
     pixels will be squared
          $ astarithmetic image.fits 2 pow

‘sqrt’
     The square root of the first operand, so "‘5 sqrt’" is equivalent
     to $\sqrt{5}$.  Since the square root is only defined for positive
     values, any negative-valued pixel will become NaN (blank).  The
     output will have a floating point type, but its precision is
     determined from the input: if the input is a 64-bit floating point,
     the output will also be 64-bit.  Otherwise, the output will be
     32-bit floating point (see *note Numeric data types:: for the
     respective precision).  Therefore if you require 64-bit precision
     in estimating the square root, convert the input to 64-bit floating
     point first, for example, with ‘5 float64 sqrt’.  For example, each
     pixel of the output of the command below will be the square root of
     that pixel in the input.
          $ astarithmetic image.fits sqrt

     If you just want to scale an image with negative values using this
     operator (for better visual inspection, and the actual values do
     not matter for you), you can subtract the image from its minimum
     value, then take its square root:

          $ astarithmetic image.fits image.fits minvalue - sqrt -g1

     Alternatively, to avoid reading the image into memory two times,
     you can use the ‘set-’ operator to read it into the variable ‘i’
     and use ‘i’ two times to speed up the operation (described below):

          $ astarithmetic image.fits set-i i i minvalue - sqrt

‘log’
     Natural logarithm of first operand, so "‘4 log’" is equivalent to
     $ln(4)$.  Negative pixels will become NaN, and the output type is
     determined from the input, see the explanation under ‘sqrt’ for
     more on these features.  For example, the command below will take
     the natural logarithm of every pixel in the input.

          $ astarithmetic image.fits log --output=log.fits

‘log10’
     Base-10 logarithm of first popped operand, so "‘4 log10’" is
     equivalent to $log_{10}(4)$.  Negative pixels will become NaN, and
     the output type is determined from the input, see the explanation
     under ‘sqrt’ for more on these features.  For example, the command
     below will take the base-10 logarithm of every pixel in the input.

          $ astarithmetic image.fits log10


File: gnuastro.info,  Node: Trigonometric and hyperbolic operators,  Next: Constants,  Prev: Basic mathematical operators,  Up: Arithmetic operators

6.2.4.2 Trigonometric and hyperbolic operators
..............................................

All the trigonometric and hyperbolic functions are described here.  One
good thing with these operators is that they take inputs and outputs in
degrees (which we usually need as input or output), not radians (like
most other programs/libraries).

‘sin’
‘cos’
‘tan’
     Basic trigonometric functions.  They take one operand, in units of
     degrees.

‘asin’
‘acos’
‘atan’
     Inverse trigonometric functions.  They take one operand and the
     returned values are in units of degrees.

‘atan2’
     Inverse tangent (output in units of degrees) that uses the signs of
     the input coordinates to distinguish between the quadrants.  This
     operator therefore needs two operands: the first popped operand is
     assumed to be the X axis position of the point, and the second
     popped operand is its Y axis coordinate.

     For example, see the commands below.  To be more clear, we are
     using Table's *note Column arithmetic:: which uses exactly the same
     internal library function as the Arithmetic program for images.  We
     are showing the results for four points in the four quadrants of
     the 2D space (if you want to try running them, you do not need to
     type/copy the parts after <#>).  The first point (2,2) is in the
     first quadrant, therefore the returned angle is 45 degrees.  But
     the second, third and fourth points are in the quadrants of the
     same order, and the returned angles reflect the quadrant.

          $ echo " 2  2" | asttable -c'arith $2 $1 atan2'   # -->   45
          $ echo " 2 -2" | asttable -c'arith $2 $1 atan2'   # -->  -45
          $ echo "-2 -2" | asttable -c'arith $2 $1 atan2'   # --> -135
          $ echo "-2  2" | asttable -c'arith $2 $1 atan2'   # -->  135

     However, if you simply use the classic arc-tangent operator
     (‘atan’) for the same points, the result will only be in two
     quadrants as you see below:

          $ echo " 2  2" | asttable -c'arith $2 $1 / atan'  # -->   45
          $ echo " 2 -2" | asttable -c'arith $2 $1 / atan'  # -->  -45
          $ echo "-2 -2" | asttable -c'arith $2 $1 / atan'  # -->   45
          $ echo "-2  2" | asttable -c'arith $2 $1 / atan'  # -->  -45

‘sinh’
‘cosh’
‘tanh’
     Hyperbolic sine, cosine, and tangent.  These operators take a
     single operand.

‘asinh’
‘acosh’
‘atanh’
     Inverse Hyperbolic sine, cosine, and tangent.  These operators take
     a single operand.


File: gnuastro.info,  Node: Constants,  Next: Coordinate conversion operators,  Prev: Trigonometric and hyperbolic operators,  Up: Arithmetic operators

6.2.4.3 Constants
.................

During your analysis it is often necessary to have certain constants
like the number $\pi$.  The "operators" in this section do not actually
take any operand, they just replace the desired constant into the stack.
So in effect, these are actually operands.  But since their value is not
inserted by the user, we have placed them in the list of operators.

‘e’
     Euler’s number, or the base of the natural logarithm (no units).
     See Wikipedia
     (https://en.wikipedia.org/wiki/E_(mathematical_constant)).

‘pi’
     Ratio of circle’s circumference to its diameter (no units).  See
     Wikipedia (https://en.wikipedia.org/wiki/Pi).

‘c’
     The speed of light in vacuum, in units of $m/s$.  see Wikipedia
     (https://en.wikipedia.org/wiki/Speed_of_light).

‘G’
     The gravitational constant, in units of $m^3/kg/s^2$.  See
     Wikipedia (https://en.wikipedia.org/wiki/Gravitational_constant).

‘h’
     Plank's constant, in units of $J/Hz$ or $kg\times m^2/s$.  See
     Wikipedia (https://en.wikipedia.org/wiki/Planck_constant).

‘au’
     Astronomical Unit, in units of meters.  See Wikipedia
     (https://en.wikipedia.org/wiki/Astronomical_unit).

‘ly’
     Distance covered by light in vacuum in one year, in units of
     meters.  See Wikipedia (https://en.wikipedia.org/wiki/Light-year).

‘avogadro’
     Avogadro's constant, in units of $1/mol$.  See Wikipedia
     (https://en.wikipedia.org/wiki/Avogadro_constant).

‘fine-structure’
     The fine-structure constant (no units).  See Wikipedia
     (https://en.wikipedia.org/wiki/Fine-structure_constant).


File: gnuastro.info,  Node: Coordinate conversion operators,  Next: Unit conversion operators,  Prev: Constants,  Up: Arithmetic operators

6.2.4.4 Coordinate conversion operators
.......................................

Different celestial coordinate systems are useful for different
scenarios.  For example, assume you have the RA and Dec of large sample
of galaxies that you plan to study the halos of galaxies from.  For such
studies, you prefer to stay as far away as possible from the Galactic
plane, because the density of stars and interstellar filaments (cirrus)
significantly increases as you get close to the Milky way's disk.  But
the Equatorial coordinate system
(https://en.wikipedia.org/wiki/Equatorial_coordinate_system) which
defines the RA and Dec and is based on Earth's equator; and does not
show the position of your objects in relation to the galactic disk.

   The best way forward in the example above is to convert your RA and
Dec table into the Galactic coordinate system
(https://en.wikipedia.org/wiki/Galactic_coordinate_system); and select
those with a large (positive or negative) Galactic latitude.
Alternatively, if you observe a bright point on a galaxy and want to
confirm if it was actually a super-nova and not a moving asteroid, a
first step is to convert your RA and Dec to the Ecliptic coordinate
system (https://en.wikipedia.org/wiki/Ecliptic_coordinate_system) and
confirm if you are sufficiently distant from the ecliptic (plane of the
Solar System; where fast moving objects are most common).

   The operators described in this section are precisely for the purpose
above: to convert various celestial coordinate systems that are
supported within Gnuastro into each other.  For example, if you want to
convert the RA and Dec equatorial (at the Julian year 2000 equinox)
coordinates (within the ‘RA’ and ‘DEC’ columns) of ‘points.fits’ into
Galactic longitude and latitude, you can use the command below (the
column metadata are not mandatory, but to avoid later confusion, it is
always good to have them in your output.

     $ asttable points.fits -c'arith RA DEC eq-j2000-to-galactic' \
                --colmetadata=1,GLON,deg,"Galactic longitude" \
                --colmetadata=2,GLAT,deg,"Galactic latitude" \
                --output=points-gal.fits

   One important thing to consider is that the equatorial and ecliptic
coordinates are not static: they include the dynamics of Earth in the
solar system: in particular, the reference point on the equator moves
over decades.  Therefore these two (equatorial and ecliptic) coordinate
systems are defined within epochs: the 1950 epoch is defined by
Besselian years
(https://en.wikipedia.org/wiki/Epoch_(astronomy)#Besselian_years), while
the 2000 epoch is defined in Julian years
(https://en.wikipedia.org/wiki/Epoch_(astronomy)#Julian_years_and_J2000).
So when dealing with these coordinates one of the '‘-b1950’' or
'‘-j2000’' suffixes are necessary (for example ‘eq-j2000’ or
‘ec-b1950’).

   The Galactic or Supergalactic coordinates are not defined based on
the Earth's dynamics; therefore they do not have any epoch associated
with them.  Extra-galactic studies do not depend on the dynamics of the
earth, but the equatorial coordinate system is the most dominant in that
field.  Therefore in its 23rd General Assembly, the International
Astronomical Union approved the International Celestial Reference System
(https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations)
or ICRS based on quasars (which are static within our observational
limitations)viewed through long baseline radio interferometry (the most
accurate method of observation that we currently have).  ICRS is
designed to be within the errors of the Equatorial J2000 coordinate
system, so they are currently very similar; but ICRS has much better
accuracy.  We will be adding ICRS in the operators below soon.

   *Floating point errors:* The operation to convert between the
coordinate systems involves many sines, cosines (and their inverse).
Therefore, floating point errors (due to the limited precision of the
definition of floating points in bits) can cause small offsets.  For
example see the code below were we convert equatorial to galactic and
back, then compare the input and output (which is in the 5th and 6th
decimal of a degree; or about 0.2 or 0.01 arcseconds).

     $ sys1=eq-j2000
     $ sys2=galactic
     $ echo "10.2345689 45.6789012" \
            | asttable -Afixed -B8 \
                       -c'arith $1 $2 '$sys1'-to-'$sys2' \
                                '$sys2'-to-'$sys1' set-lat set-lng \
                                lng $1 - lat $2 -'
     0.00000363   -0.00007725

   If you set ‘sys2=ec-j2000’ or ‘sys2=supergalactic’, it will be zero
until the full set of 8 decimals that are printed here (the displayed
precision can be changed with the value of the ‘-B’ option above).  It
is therefore useful to have your original coordinates (in the same table
for example) and not do too many conversions on conversions (to
propagate this problem).

‘eq-b1950-to-eq-j2000’
‘eq-b1950-to-ec-b1950’
‘eq-b1950-to-ec-j2000’
‘eq-b1950-to-galactic’
‘eq-b1950-to-supergalactic’
     Convert Equatorial (B1950 equinox) coordinates into the respective
     coordinate system within each operator's name.

‘eq-j2000-to-eq-b1950’
‘eq-j2000-to-ec-b1950’
‘eq-j2000-to-ec-j2000’
‘eq-j2000-to-galactic’
‘eq-j2000-to-supergalactic’
     Convert Equatorial (J2000 equinox) coordinates into the respective
     coordinate system within each operator's name.

‘ec-b1950-to-eq-b1950’
‘ec-b1950-to-eq-j2000’
‘ec-b1950-to-ec-j2000’
‘ec-b1950-to-galactic’
‘ec-b1950-to-supergalactic’
     Convert Ecliptic (B1950 equinox) coordinates into the respective
     coordinate system within each operator's name.

‘ec-j2000-to-eq-b1950’
‘ec-j2000-to-eq-j2000’
‘ec-j2000-to-ec-b1950’
‘ec-j2000-to-galactic’
‘ec-j2000-to-supergalactic’
     Convert Ecliptic (J2000 equinox) coordinates into the respective
     coordinate system within each operator's name.

‘galactic-to-eq-b1950’
‘galactic-to-eq-j2000’
‘galactic-to-ec-b1950’
‘galactic-to-ec-j2000’
‘galactic-to-supergalactic’
     Convert Galactic coordinates into the respective coordinate system
     within each operator's name.

‘supergalactic-to-eq-b1950’
‘supergalactic-to-eq-j2000’
‘supergalactic-to-ec-b1950’
‘supergalactic-to-ec-j2000’
‘supergalactic-to-galactic’
     Convert Supergalactic coordinates into the respective coordinate
     system within each operator's name.


File: gnuastro.info,  Node: Unit conversion operators,  Next: Statistical operators,  Prev: Coordinate conversion operators,  Up: Arithmetic operators

6.2.4.5 Unit conversion operators
.................................

It often happens that you have data in one unit (for example, counts on
your CCD), but would like to convert it into another (for example,
magnitudes, to measure the brightness of a galaxy).  While the equations
for the unit conversions can be easily found on the internet, the
operators in this section are designed to simplify the process and let
you do it easily and fast without having to remember constants and
relations.

‘counts-to-mag’
     Convert counts (usually CCD outputs) to magnitudes using the given
     zero point.  The zero point is the first popped operand and the
     count image or value is the second popped operand.

     For example, assume you have measured the standard deviation of the
     noise in an image to be ‘0.1’ counts, and the image's zero point is
     ‘22.5’ and you want to measure the _per-pixel_ surface brightness
     limit of the dataset(1).  To apply this operator on an image,
     simply replace ‘0.1’ with the image name, as described below.

          $ astarithmetic 0.1 22.5 counts-to-mag --quiet

     Of course, you can also convert every pixel in an image (or table
     column in Table's *note Column arithmetic::) with this operator if
     you replace the second popped operand with an image/column name.
     For an example of applying this operator on an image, see the
     description of surface brightness in *note Brightness flux
     magnitude::, where we will convert an image's pixel values to
     surface brightness.

‘mag-to-counts’
     Convert magnitudes to counts (usually CCD outputs) using the given
     zero point.  The zero point is the first popped operand and the
     magnitude value is the second.  For example, if an object has a
     magnitude of 20, you can estimate the counts corresponding to it
     (when the image has a zero point of 24.8) with this command: Note
     that because the output is a single number, we are using ‘--quiet’
     to avoid printing extra information.

          $ astarithmetic 20 24.8 mag-to-counts --quiet

‘mag-to-luminosity’
     Convert the given apparent magnitude to luminosity (in units of the
     luminosity of a reference object).  It takes the following three
     operands (in the order written on the command-line; see example
     below):

       1. The measured apparent magnitude that is corrected for the ISM
          absorption/extinction.  You can find the ISM
          absorption/extinction of a certain position on the sky using
          Gnuastro's Query program, which gives you direct access to the
          NED Extinction Calculator as described in *note Available
          databases::.

       2. Absolute magnitude of the reference object in the same filter
          and magnitude system (for example Vega or AB) that the
          measured apparent magnitude was derived from.  The reference
          object is conventionally the Sun.  The Sun's absolute
          magnitude in various commonly used filters and magnitude
          systems from table 3 of Willmer 2018
          (https://arxiv.org/abs/1804.07788).  For example the Sun's
          absolute magnitude in the SDSS u, g, r, i and z filters (with
          the AB magnitude system) is respectively 6.39, 5.11, 4.65,
          4.53 and 4.50.
       3. Difference of apparent ($m$) and absolute ($M$) magnitudes:
          $m-M$.  At small distances or when the input is bolometric
          (across all wavelengths), the distance modulus can be used for
          this.  See ‘--distancemodulus’ in *note CosmicCalculator basic
          cosmology calculations:: for more details and why
          ‘--absmagconv’ is preferred to the distance modulus in the
          absense of SED-based methods.

     For example, let's assume the apparent AB magnitude of a galaxy
     (after correcting Galactic extinction) is 20 in the SDSS g filter,
     and that its redshift is 0.01 and we need its luminosity (in units
     of solar luminosity) in the same filter.  To do this, we need the
     apparent-absolute magnitude difference (using CosmicCalculator), as
     well as the Sun's absolute magnitude, from Willmer 2018
     (https://arxiv.org/abs/1804.07788) in this filter (which is 5.11):

          $ conv=$(astcosmiccal --absmagconv --redshift=0.01)
          $ astarithmetic 20 5.11 $conv mag-to-luminosity

‘luminosity-to-mag’
     Convert luminostiy (in units of solar luminosity) to apparent
     magnitude using the distance modulus.  This is the inverse of
     ‘mag-to-luminosity’, see its description for the details and usage
     example.

‘counts-to-sb’
     Convert counts to surface brightness using the zero point and area
     (in units of arcsec$^2$).  The first popped operand is the area (in
     arcsec$^2$), the second popped operand is the zero point and the
     third are the count values.  Estimating the surface brightness
     involves taking the logarithm.  Therefore this operator will
     produce NaN for counts with a negative value.

     For example, with the commands below, we read the zero point from
     the image headers (assuming it is in the ‘ZPOINT’ keyword), we
     calculate the pixel area from the image itself, and we call this
     operator to convert the image pixels (in counts) to surface
     brightness (mag/arcsec$^2$).

          $ zeropoint=$(astfits image.fits --keyvalue=ZPOINT -q)
          $ pixarea=$(astfits image.fits --pixelareaarcsec2)
          $ astarithmetic image.fits $zeropoint $pixarea counts-to-sb \
                          --output=image-sb.fits

     For more on the definition of surface brightness see *note
     Brightness flux magnitude::, and for a fully tutorial on optimal
     usage of this, see *note FITS images in a publication::.

‘sb-to-counts’
     Convert surface brightness using the zero point and area (in units
     of arcsec$^2$) to counts.  The first popped operand is the area (in
     arcsec$^2$), the second popped operand is the zero point and the
     third are the surface brightness values.  See the description of
     ‘counts-to-sb’ for more.

‘mag-to-sb’
     Convert magnitudes to surface brightness over a certain area (in
     units of arcsec$^2$).  The first popped operand is the area and the
     second is the magnitude.  For example, let's assume you have a
     table with the two columns of magnitude (called ‘MAG’) and area
     (called ‘AREAARCSEC2’).  In the command below, we will use *note
     Column arithmetic:: to return the surface brightness.
          $ asttable table.fits -c'arith MAG AREAARCSEC2 mag-to-sb'

‘sb-to-mag’
     Convert surface brightness to magnitudes over a certain area (in
     units of arcsec$^2$).  The first popped operand is the area and the
     second is the magnitude.  See the description of ‘mag-to-sb’ for
     more.

‘counts-to-jy’
     Convert counts (usually CCD outputs) to Janskys through an
     AB-magnitude based zero point.  The top-popped operand is assumed
     to be the AB-magnitude zero point and the second-popped operand is
     assumed to be a dataset in units of counts (an image in Arithmetic,
     and a column in Table's *note Column arithmetic::).  For the full
     equation and basic definitions, see *note Brightness flux
     magnitude::.

     For example, SDSS images are calibrated in units of nanomaggies,
     with a fixed zero point magnitude of 22.5.  Therefore you can
     convert the units of SDSS image pixels to Janskys with the command
     below:

          $ astarithmetic sdss-image.fits 22.5 counts-to-jy

‘jy-to-counts’
     Convert Janskys to counts (usually CCD outputs) through an
     AB-magnitude based zero point.  This is the inverse operation of
     the ‘counts-to-jy’, see there for usage example.

‘zeropoint-change’
     Change the zero point of the input dataset to a new zero point
     (output will always be 64-bit floating point).  For example, with
     the command below, we are changing the zero point of an image from
     26.892 to 27.0:

          $ astarithmetic image.fits 26.892 27.0 zeropoint-change

     Note that the input or output zero points can also be an image
     (with the same size as the input image).  This is very useful to
     correct situations where the zero point can change over the image
     (happens in single exposures) and you want to unify it across the
     whole image (to build a deep coadded image).

‘counts-to-nanomaggy’
     Convert counts to Nanomaggy (with fixed zero point of 22.5, used as
     the pixel units of many surveys like SDSS). For example if your
     image has a zero point of 24.93, you can convert it to Nanomaggies
     with the command below:

          $ astarithmetic image.fits 24.93 counts-to-nanomaggy

     This is just a wrapper over the ‘zeropoint-change’ operator where
     the output zero point is 22.5.

‘nanomaggy-to-counts’
     Convert Nanomaggy to counts.  Nanomaggy is defined to have a fixed
     zero point of 22.5 and is the pixel units of many surveys like
     SDSS. For example if you would like to convert an image in units of
     Nanomaggy (for example from SDSS) to the counts of a camera with a
     zero point of 25.92, you can use the command below:

          $ astarithmetic image.fits 25.92 nanomaggy-to-counts

     This is just a wrapper over the ‘zeropoint-change’ operator where
     the input zero point is 22.5.

‘mag-to-jy’
     Convert AB magnitudes to Janskys, see *note Brightness flux
     magnitude::.

‘jy-to-mag’
     Convert Janskys to AB magnitude, see *note Brightness flux
     magnitude::.

‘jy-to-wavelength-flux-density’
     Convert Janskys to wavelength flux density (in units of
     $erg/cm^2/s/\AA$) at a certain wavelength (given in units of
     Angstroms).  Recall that Jansky is also a unit of spectral flux
     density, but it is in units of frequency ($erg/cm^2/s/Hz$).

     For example at the wavelength of 5556, Vega's frequency flux
     density is 3500 Janskys.  To convert this to wavelength flux
     density, you can use this command:

          $ astarithmetic 3.5e3 5556 jy-to-wavelength-flux-density

     If your input values are in units of counts or magnitudes, you can
     use the ‘mag-to-jy’ or ‘count-to-jy’ operators that are described
     above.  For example, if you want the wavelength flux density of a
     source with magnitude 21 in the J-PAS J0660 filter (centered at
     approximately $6600\AA$), you can use this command:

          $ astarithmetic 21 mag-to-jy 6600 \
                          jy-to-wavelength-flux-density

     The conversion is done based on this derivation: the speed of light
     ($c$) can be written as:

    $$c=2.99792458\times10^8 m/s = 2.99792458\times10^{18} \AA Hz$$

     The speed of light connects the wavelength ($\lambda$) and
     frequency ($\nu$) of photons: $c=\nu\lambda$.  Therefore,
     $\nu=c/\lambda$ and taking the derivative:
     $d\nu=(c/\lambda^2)d\lambda$ or $d\nu/d\lambda=c/\lambda^2$.
     Inserting physical values and units:

$$d\nu/d\lambda = \frac{2.99792458\times10^{18} \AA Hz}{\lambda^2\AA^2} = \frac{2.99792458\times10^{18}}{\lambda^2} Hz/A$$

     Recall that to convert a function of $\nu$ into a function of
     $\lambda$, where $\nu$ and $\lambda$ are also related to each
     other, we have the following equation: $f(\lambda) =
     \frac{d\nu}{d\lambda} f(\nu)$.  Here, $f(\nu)$ is the value in
     Janskys ($J\times10^{-23}erg/s/cm^2/Hz$; see *note Brightness flux
     magnitude::).  Replacing $d\nu/d\lambda$ from above we get:

$$F_\lambda=\frac{2.99792458\times10^8}{\lambda^2}Hz/A \times J\times10^{-23}erg/s/cm^2/Hz$$

$$F_\lambda=J\times\frac{2.99792458\times10^{-5}}{\lambda^2} erg/s/cm^2/\AA$$

‘wavelength-flux-density-to-jy’
     Convert wavelength flux density (in units of $erg/cm^2/s/\AA$) to
     Janskys at a certain wavelength (given in units of Angstroms).  For
     details and usage examples, see the description of
     ‘jy-to-wavelength-flux-density’ (the inverse of this function).

‘sblim-diff’
     Calculate the difference in surface brightness limit after using a
     different exposed radius (of a different telescope) or a different
     exposure time.  The filter used, the multiple of sigma and the
     extrapolated area of the derived surface brightness limit will be
     the same as your reference surface brighntess limit.  For more on
     the surface brightness limit, the exposed radius, the extrapolated
     area and the equation behind this operator, see *note Surface
     brightness limit of image::.

     For example, we know that in the rSDSS filter, the SDSS survey
     reaches a $3\sigma$ surface brightness limit of 26.5 mag/arcsec$^2$
     over 100 arcsec$^2$.  Let's use this to find the expected surface
     brighntess limit of the LSST survey in the same filter (after its
     observations have been completed).  To do this, we need to know the
     exposed radius and exposure time of SDSS and LSST to feed into the
     command below:
        • The SDSS telescope(2) has an exposed radius of
          $\sqrt{1.25^2-0.54^2}=1.13$ meters and it had a 1 minute (60
          seconds) exposure time.
        • The Vera C. Rubin telescope(3) (used for LSST) has an exposed
          radius of $\sqrt{4.2^2-2.5^2}=3.37$ meters and at its end, it
          will have 8 hours (28800 seconds) of total exposure on all
          observed regions(4) (in all filters).  Because it has 6
          filters, in the rSDSS filter it will have a total exposure
          time of about $28800/6=4800$ seconds in its observed regions.

          $ astarithmetic 3.37 1.13 / 4800 60 / sblim-diff 26.5 +
          3.00652410811219e+01
          Arithmetic finished in 0.012210 seconds

     To have an easy to read number, you can pipe it to the Table
     program and use the ‘-Y’ option.  But the time-reporting feature of
     Arithmetic will be problematic, so you should also pass ‘--quiet’
     (or ‘-q’) to Arithmetic.

          $ astarithmetic 3.37 1.13 / 4800 60 / sblim-diff 26.5 + -q \
                          | asttable -Y
          30.065241

     Therefore, juduging purely based on exposed area and time (too
     idealistic!), the LSST will reach a $3\sigma$ surface brighntess
     limit of 30.07 over 100 arcsec$^2$.

     A simple note for the second command: this trick is only useful for
     operators like this that are sure to produce numbers that are
     larger than $10^{-2}$.  You will loose precision for values smaller
     than that and using the default scientific notation is best.

‘au-to-pc’
     Convert Astronomical Units (AUs) to Parsecs (PCs).  This operator
     takes a single argument which is interpreted to be the input AUs.
     The conversion is based on the definition of Parsecs: $1 \rm{PC} =
     1/tan(1^{\prime\prime}) \rm{AU}$, where $1^{\prime\prime}$ is one
     arcseconds.  In other words, $1 (\rm{PC}) = 648000/\pi (\rm{AU})$.
     For example, if we take Pluto's average distance to the Sun to be
     40 AUs, we can obtain its distance in Parsecs using this command:

          echo 40 | asttable -c'arith $1 au-to-pc'

‘pc-to-au’
     Convert Parsecs (PCs) to Astronomical Units (AUs).  This operator
     takes a single argument which is interpreted to be the input PCs.
     For more on the conversion equation, see description of ‘au-to-pc’.
     For example, Proxima Centauri (the nearest star to the Solar
     system) is 1.3020 Parsecs from the Sun, we can calculate this
     distance in units of AUs with the command below:

          echo 1.3020 | asttable -c'arith $1 pc-to-au'

‘ly-to-pc’
     Convert Light-years (LY) to Parsecs (PCs).  This operator takes a
     single argument which is interpreted to be the input LYs.  The
     conversion is done from IAU's definition of the light-year
     (9460730472580800 m $\approx$ 63241.077 AU = 0.306601 PC, for the
     conversion of AU to PC, see the description of ‘au-to-pc’).

     For example, the distance of Andromeda galaxy to our galaxy is 2.5
     million light-years, so its distance in kilo-Parsecs can be
     calculated with the command below (note that we want the output in
     kilo-parsecs, so we are dividing the output of this operator by
     1000):

          echo 2.5e6 | asttable -c'arith $1 ly-to-pc 1000 /'

‘pc-to-ly’
     Convert Parsecs (PCs) to Light-years (LY). This operator takes a
     single argument which is interpreted to be the input PCs.  For the
     conversion and an example of the inverse of this operator, see the
     description of ‘ly-to-pc’.

‘ly-to-au’
     Convert Light-years (LY) to Astronomical Units (AUs).  This
     operator takes a single argument which is interpreted to be the
     input LYs.  For the conversion and a similar example, see the
     description of ‘ly-to-pc’.

‘au-to-ly’
     Convert Astronomical Units (AUs) to Light-years (LY). This operator
     takes a single argument which is interpreted to be the input AUs.
     For the conversion and a similar example, see the description of
     ‘ly-to-pc’.

   ---------- Footnotes ----------

   (1) The _per-pixel_ surface brightness limit is the magnitude of the
noise standard deviation.  For more on surface brightness see *note
Brightness flux magnitude::.  In the example command, because the output
is a single number, we are using ‘--quiet’ to avoid printing extra
information.

   (2) 
<https://voyages.sdss.org/preflight/capturing-recording-light/sdss-telescope>

   (3) https://www.lsst.org/sites/default/files/docs/ivezic_086.02.pdf

   (4) <https://www.lsst.org/sites/default/files/docs/ivezic_086.02.pdf>


File: gnuastro.info,  Node: Statistical operators,  Next: Coadding operators,  Prev: Unit conversion operators,  Up: Arithmetic operators

6.2.4.6 Statistical operators
.............................

The operators in this section take a single dataset as input, and will
return the desired statistic as a single value.

‘minvalue’
     Minimum value in the first popped operand, so "‘a.fits minvalue’"
     will push the minimum pixel value in this image onto the stack.
     When this operator acts on a single image, the output (operand that
     is put back on the stack) will no longer be an image, but a number.
     The output of this operand is in the same type as the input.  This
     operator is mainly intended for multi-element datasets (for
     example, images or data cubes), if the popped operand is a number,
     it will just return it without any change.

     Note that when the final remaining/output operand is a single
     number, it is printed onto the standard output.  For example, with
     the command below the minimum pixel value in ‘image.fits’ will be
     printed in the terminal:
          $ astarithmetic image.fits minvalue

     However, the output above also includes a lot of extra information
     that are not relevant in this context.  If you just want the final
     number, run Arithmetic in quiet mode:
          $ astarithmetic image.fits minvalue -q

     Also see the description of ‘sqrt’ for other example usages of this
     operator.

‘maxvalue’
     Maximum value of first operand in the same type, similar to
     ‘minvalue’, see the description there for more.  For example
          $ astarithmetic image.fits maxvalue -q

‘numbervalue’
     Number of non-blank elements in first operand in the ‘uint64’ type
     (since it is always a positive integer, see *note Numeric data
     types::).  Its usage is similar to ‘minvalue’, for example
          $ astarithmetic image.fits numbervalue -q

‘sumvalue’
     Sum of non-blank elements in first operand in the ‘float32’ type.
     Its usage is similar to ‘minvalue’, for example
          $ astarithmetic image.fits sumvalue -q

‘meanvalue’
     Mean value of non-blank elements in first operand in the ‘float32’
     type.  Its usage is similar to ‘minvalue’, for example
          $ astarithmetic image.fits meanvalue -q

‘stdvalue’
     Standard deviation of non-blank elements in first operand in the
     ‘float32’ type.  Its usage is similar to ‘minvalue’, for example
          $ astarithmetic image.fits stdvalue -q

‘medianvalue’
     Median of non-blank elements in first operand with the same type.
     Its usage is similar to ‘minvalue’, for example
          $ astarithmetic image.fits medianvalue -q

‘label-area’
     Calculate the area (number of pixels) that have a unique label and
     replace that label's pixels with its area.  The output will be an
     unsigned 32-bit integer (see *note Numeric data types::).

     The input labled image should have an integer type smaller than a
     signed 32-bit integer (so unsigned 32-bit or 64-bit integers are
     not acceptable).  For exapmle the output of ‘connected-components’
     (see *note Mathematical morphology operators::), or Gnuastro's
     *note Segment:: program.  This operator can be useful when you want
     to separate the different labels by their area.  For example let's
     assume you want to remove all labels that are smaller than 3
     pixels.  With the command below you can do that and generate new
     labels (starting from one).

          $ astarithmeitc label.fits label-area \
                          3 gt 2 connected-components

‘label-minimum’
‘label-maximum’
     Calculate the minimum/maximum value of first popped operand (values
     image) that have a unique label in the second popped operand
     (labels image).  The output will have the same type as the values
     image and the input labeled image should have the same properties
     as the ‘label-area’ operator described above.

     For example, Let's assume you want to find "thinner" objects in the
     image (to keep or delete them).  You can use the command below
     (using ‘number-neighbors’ of *note Mathematical morphology
     operators::) as the second popped operand of this operator and then
     apply a threshold of 6 for example and generate new labels for the
     remaining objects.

          astarithmetic lab.fits set-l \
                        l l 0 gt 2 number-neighbors label-maximum \
                        6 uint8 gt 2 connected-components

‘madclip-maskfilled’
‘sigclip-maskfilled’
     Mask (set to blank/NaN) all the outlying elements (defined by
     $\sigma$ or MAD clipping) in the inputs and put all the inputs back
     on the stack.  The first popped operand is the termination criteria
     of the clipping, the second popped operand is the multiple of
     $\sigma$ or MAD and the third is the number of input datasets that
     will be popped for the actual operation.  If you are not yet
     familiar with $\sigma$ or MAD clipping, it is recommended to read
     this tutorial: *note Clipping outliers::.

     When more than 95$\%$ of the area of an operand is masked, the full
     operand will be masked.  This is necessary in like this: one of
     your inputs has many outliers (for example it is much more noisy
     than the rest or its sky level has not been subtracted properly).
     Because this operator fills holes between outlying pixels, most of
     the area of the input will be masked, but the thin edges (where
     there are no "holes") will remain, causing different statistics in
     those thin edges of that input in your final coadd.  Through this
     mask coverage fraction (which is currently hard-coded(1)), we
     ensure that such thin edges do not cause artifacts in the final
     coadd.

     For example, with the second command below, we are masking the MAD
     clipped pixels of the 9 inputs (that are generated in *note
     Clipping outliers::) and writing them as separate HDUs of the
     output.  The clipping is done with 5 times the MAD and the clipping
     starts when the relative difference between subsequent MADs is
     0.01.  Finally, with the third command, we see 10 HDUs in the
     output (because the first, or 0-th, is just metadata).

          $ ls in-*.fits
          in-1.fits  in-3.fits  in-5.fits  in-7.fits  in-9.fits
          in-2.fits  in-4.fits  in-6.fits  in-8.fits

          $ astarithmetic in-*.fits 9 5 0.01 madclip-maskfilled \
                          -g1 --writeall --output=clipped.fits

          $ astfits clipped.fits --numhdus
          10

     In the Arithmetic command above, ‘--writeall’ is necessary because
     this operator puts all its inputs back on the stack of operands.
     This is because these are usually just intermediate operators.  For
     example, after masking the outliers from each input, you may want
     to coadd them into one deeper image (with the *note Coadding
     operators::.  After the coadding is done, only one operand will be
     on the stack, and ‘--writeall’ will no longer be necessary.  For
     example if you want to see how many images were used in the final
     coadd's pixels, you can use the ‘number’ operator like below:

          $ astarithmetic in-*.fits 9 5 0.01 madclip-maskfilled \
                          9 number -g1 --output=num-good.fits

‘unique’
     Remove all duplicate (and blank) elements from the first popped
     operand.  The unique elements of the dataset will be stored in a
     single-dimensional dataset.

     Recall that by default, single-dimensional datasets are stored as a
     table column in the output.  But you can use ‘--onedasimage’ or
     ‘--onedonstdout’ to respectively store them as a single-dimensional
     FITS array/image, or to print them on the standard output.

     Although you can use this operator on the floating point dataset,
     due to floating-point errors it may give non-reasonable values:
     because the tenth digit of the decimal point is also considered
     although it may be statistically meaningless, see *note Numeric
     data types::.  It is therefore better/recommended to use it on the
     integer dataset like the labeled images of *note Segment output::
     where each pixel has the integer label of the object/clump it is
     associated with.  For example, let's assume you have cropped a
     region of a larger labeled image and want to find the
     labels/objects that are within the crop.  With this operator, this
     job is trivial:
          $ astarithmetic seg-crop.fits unique

‘noblank’
     Remove all blank elements from the first popped operand.  Since the
     blank pixels are being removed, the output dataset will always be
     single-dimensional, independent of the dimensionality of the input.

     Recall that by default, single-dimensional datasets are stored as a
     table column in the output.  But you can use ‘--onedasimage’ or
     ‘--onedonstdout’ to respectively store them as a single-dimensional
     FITS array/image, or to print them on the standard output.

     For example, with the command below, the non-blank pixel values of
     ‘cropped.fits’ are printed on the command-line (the ‘--quiet’
     option is used to remove the extra information that Arithmetic
     prints as it reads the inputs, its version and its running time).

          $ astarithmetic cropped.fits noblank --onedonstdout --quiet

   ---------- Footnotes ----------

   (1) Please get in touch with us at ‘bug-gnuastro@gnu.org’ if you
notice this problem and feel the fraction needs to be lowered (or
generally to be set in each run).


File: gnuastro.info,  Node: Coadding operators,  Next: Filtering operators,  Prev: Statistical operators,  Up: Arithmetic operators

6.2.4.7 Coadding operators
..........................

The operators in this section are used when you have multiple datasets
that you would like to merge into one.  For example, you have taken ten
exposures of your scientific target, and you would like to combine them
all into one deep stacked image that is deeper.  This is commonly known
as "stacking" or "coaddition".  We use the latter in Gnuastro because
"stack" refers to the intermediate 3D data set (if we are coadding 2D
images).  However, the final product of this operation has the same
number of dimensions as each of the inputs.  The astronomical community
has traditionally used the term "coadd" to specify that the nature of
the output 2D image is different from each of the input 2D images (it is
created from them).  Furthermore, within Arithmetic, "Stack" is already
reserved for the stack of operands that the operators read from (see
*note Reverse polish notation::).

*Masking outliers (before coadding):* Outliers in one of the inputs (for
example star ghosts, satellite trails, or cosmic rays, can leave their
imprints in the final coadd.  One good way to remove them is the
‘madclip-maskfilled’ operator that can be called before the operators
here.  It is described in *note Statistical operators::; and a full
tutorial on understanding outliers and how best to remove them is
available in *note Clipping outliers::.

*Hundreds or thousands of images to coadd:* It can happen that you need
to coadd hundreds or thousands of images.  Added with the possibly long
file/directory names, this can lead to an extremely long shell command
that may cause an "Argument list too long" error in your shell.  To
avoid this, you should use Arithmetic's ‘--arguments’ option, see *note
Invoking astarithmetic::.

   When calling the coadding operators you should determine how many
operands they should take in: unlike the rest of the operators that have
a fixed number of input operands, these operators have a variable number
of input operators.  As described in the first operand below, you do
this through their first (or early) popped operands (which should be a
single integer number that is larger than one).  Below are some
important points for all the coadding operators described in this
section:

   • NaN/blank pixels will be ignored, see *note Blank pixels::.

   • The operation will be multi-threaded, greatly speeding up the
     process if you have large and numerous data to coadd.  You can
     disable multi-threaded operations with the ‘--numthreads=1’ option
     (see *note Multi-threaded operations::).

‘min’
‘max’
‘sum’
‘std’
‘mad’
‘mean’
‘median’
‘number’
     For each pixel, calculate the respective statistic from in all
     given datasets.  For the ‘min’ and ‘max’ operators, the output will
     have the same type as the input.  For the ‘number’ operator, the
     output will have an unsigned 32-bit integer type and the rest will
     be 32-bit floating point.

     The first popped operand to this operator must be a positive
     integer number which specifies how many further operands should be
     popped from the stack.  All the subsequently popped operands must
     have the same type and size.  This operator (and all the
     variable-operand operators similar to it that are discussed below)
     will work in multi-threaded mode unless Arithmetic is called with
     the ‘--numthreads=1’ option, see *note Multi-threaded operations::.

     For example, the following command will produce an image with the
     same size and type as the three inputs, but each output pixel value
     will be the minimum of the same pixel's values in all three input
     images.

          $ astarithmetic a.fits b.fits c.fits 3 min --output=min.fits

     Regarding the ‘number’ operator: some datasets may have blank
     values (which are also ignored in all similar operators like ‘min’,
     ‘sum’, ‘mean’ or ‘median’).  Hence, the final pixel values of this
     operator will not, in general, be equal to the number of inputs.
     This operator is therefore mostly called in parallel with those
     operators to know the "weight" of each pixel (in case you want to
     only keep pixels that had the full exposure for example).

‘quantile’
     For each pixel, find the quantile from all given datasets.  The
     output will have the same numeric data type and size as the input
     datasets.  Besides the input datasets, the quantile operator also
     needs a single parameter (the requested quantile).  The parameter
     should be the first popped operand, with a value between (and
     including) 0 and 1.  The second popped operand must be the number
     of datasets to use.

     In the example below, the first-popped operand (‘0.7’) is the
     quantile, the second-popped operand (‘3’) is the number of datasets
     to pop.

          astarithmetic a.fits b.fits c.fits 3 0.7 quantile

‘sigclip-mad’
‘sigclip-std’
‘sigclip-mean’
‘sigclip-median’
     Return the respective statistic after $\sigma$-clipping the values
     of the same pixel of all the input operands.  The respective
     statistic will be stored in a 32-bit floating point number.  The
     number of inputs used to make the desired measurement for each
     pixel is also returned as a second output operand; see below for
     more on how to deal with the second output operand.

     For a complete tutorial on clipping outliers when coadding images
     see *note Clipping outliers:: (if you haven't read it yet, we
     encourage you to read through it before continuing).  In
     particular, the most robust solution is to first use
     ‘madclip-maskfilled’ (described in *note Statistical operators::),
     then use any of these.

     This operator is very similar to ‘min’, with the exception that it
     expects two extra operands (parameters for MAD-clipping) before the
     total number of inputs.  The first popped operand is the
     termination criteria and the second is the multiple of the median
     absolute deviation.

     For example, in the command below, the first popped operand of
     ‘sigclip-mean’ (‘0.1’) is the $\sigma$-clipping termination
     criteria.  If the termination criteria is larger than, or equal to
     1, it is interpreted as the total number of clips.  But if it is
     between 0 and 1, then it is the tolerance level on the change in
     the median absolute deviation (see *note Sigma clipping::).  The
     second popped operand (‘4’) is the multiple of sigma (STD) to use.
     The third popped operand (‘3’) is number of datasets that should be
     coadded (similar to the first popped operand to ‘min’).  Two other
     side-notes should be mentioned here:
        • As mentioned above, before this operator, we are masking the
          filled MAD-clipped elements with ‘madclip-maskfilled’.  As
          described in *note Clipping outliers::, this is very important
          for removing the types of outliers that we have in
          astronomical imaging.
        • We are using ‘--writeall’ because this operator places two
          operands on the coadd: your desired statistics, and the number
          of inputs that were used in it (after clipping).

          $ astarithmetic a.fits b.fits c.fits -g1 --writeall \
                          3 5 0.01 madclip-maskfilled \
                          3 4 0.1  sigclip-mean

     The numbers image has the smallest unsigned integer type that fits
     the total number of your input datasets (see *note Numeric data
     types::).  For example if you have less than 255 input operands
     (not pixels!), then it will have an unsigned 8-bit integer type, if
     you have 1000 input operands (or any number less than 65534
     inputs), it will be an unsigned 16-bit integer.  Recall that when
     you have many input files to coadd, it may be necessary to write
     the arguments into a text file and use ‘--arguments’ (see *note
     Invoking astarithmetic::).

     The numbers image is included by default because it is usually
     important in clipping based coadds (where the number of inputs used
     in the calculation of each pixel can be different from another
     pixel, and this affects the final output noise).  In case you are
     not interested in the numbers image, you should first ‘swap’ the
     two output operands, then ‘free’ the top operand like below.

          $ astarithmetic a.fits b.fits c.fits -g1 \
                          3 5 0.01 madclip-maskfilled \
                          3 4 0.1  sigclip-mean swap free \
                          --output=single-hdu.fits

     In case you just want the numbers image, you can use
     ‘sigclip-median’ (which is always calculated as part of the
     clipping process: no extra overhead), and ‘free’ the top operand
     (without the ‘swap’: the median-coadd image), leaving only the
     numbers image:

          $ astarithmetic a.fits b.fits c.fits -g1 \
                          3 5 0.01 madclip-maskfilled \
                          3 4 0.1  sigclip-median free \
                          --output=single-hdu-only-numbers.fits

‘madclip-mad’
‘madclip-std’
‘madclip-mean’
‘madclip-median’
     Similar to the ‘sigclip-*’ operators, but using Median Absolute
     Deviation (MAD) as the measure of spread.  See *note MAD clipping::
     and more generally *note Clipping outliers::.  See the description
     of ‘sigclip-*’ for usage details; just remember that o1 MAD is
     equivalent to $0.67449\sigma$.

‘sigclip-all’
‘madclip-all’
     Similar to the ‘sigclip-*’ or ‘madclip-*’ operators, but instead of
     returning only a single measured statistic (mean, STD, median or
     MAD) all of them are returned, and in the following order: mean,
     STD, median, MAD, number.  This will greatly speed up your analysis
     pipelines when more than one statistic is required (for example you
     need both the mean and the standard deviation).  Just note that as
     described under the sigma-clip operators, ‘--writeall’ will be
     necessary if you want them all to be put in the output.

     For example, let's assume you want the mean and standard deviation
     after a mad-clipping like the examples before (not just one of
     them).  In that case, the cleanest way is to use a combination of
     the ‘set-’ and ‘free’ operators like the example below.  The first
     (‘set-’ pops the top operand from the stack and stores it in a
     named variable (that you can use any number of times later), while
     the second (‘free’) simply pops and deletes the top operand on the
     stack.  For more on these operators, see *note Operand storage in
     memory or a file::.  Finally, to help in later usage of the HDUs,
     we use the ‘--metadata’ operator to add a name to each output HDU
     in the same Arithmetic call.

          $ astarithmetic a.fits b.fits c.fits -g1 --writeall \
                          3 7.5 0.2 madclip-all set-m set-s \
                          free free free s m --metaname=MEAN,STD \
                          --output=out.fits

     Different combinations of the various statistics are commonly
     necessary for differnet purposes.  The rason we haven't specified
     an option for each combination is because that would dramatically
     increase the operator-count and would confuse the developers and
     users.  On the other hand, the process of maskfilling or even just
     clipping is much more computationally expensive on large datasets
     than these measurements combined.  We therefore recommend the
     process above to extract the statistics you want.


File: gnuastro.info,  Node: Filtering operators,  Next: Pooling operators,  Prev: Coadding operators,  Up: Arithmetic operators

6.2.4.8 Filtering (smoothing) operators
.......................................

Image filtering is commonly used for smoothing: every pixel value in the
output image is created by applying a certain statistic to the pixels in
its vicinity.

‘filter-mean’
     Apply mean filtering (or moving average
     (https://en.wikipedia.org/wiki/Moving_average)) on the input
     dataset.  During mean filtering, each pixel (data element) is
     replaced by the mean value of all its surrounding pixels (excluding
     blank values).  The number of surrounding pixels in each dimension
     (to calculate the mean) is determined through the earlier operands
     that have been pushed onto the stack prior to the input dataset.
     The number of necessary operands is determined by the dimensions of
     the input dataset (first popped operand).  The order of the
     dimensions on the command-line is the order in FITS format.  Here
     is one example:

          $ astarithmetic 5 4 image.fits filter-mean

     In this example, each pixel is replaced by the mean of a 5 by 4 box
     around it.  The box is 5 pixels along the first FITS dimension
     (horizontal when viewed in ds9) and 4 pixels along the second FITS
     dimension (vertical).

     Each pixel will be placed in the center of the box that the mean is
     calculated on.  If the given width along a dimension is even, then
     the center is assumed to be between the pixels (not in the center
     of a pixel).  When the pixel is close to the edge, the pixels of
     the box that fall outside the image are ignored.  Therefore, on the
     edge, less points will be used in calculating the mean.

     The final effect of mean filtering is to smooth the input image, it
     is essentially a convolution with a kernel that has identical
     values for all its pixels (is flat), see *note Convolution
     process::.

     Note that blank pixels will also be affected by this operator: if
     there are any non-blank elements in the box surrounding a blank
     pixel, in the filtered image, it will have the mean of the
     non-blank elements, therefore it will not be blank any more.  If
     blank elements are important for your analysis, you can use the
     ‘isblank’ operator with the ‘where’ operator to set them back to
     blank after filtering.

     For example in the command below, we are first filtering the image,
     then setting its original blank elements back to blank in the
     output of filtering (all within one Arithmetic command).  Note how
     we are using the ‘set-’ operator to give names to the temporary
     outputs of steps and simplify the code (see *note Operand storage
     in memory or a file::).

          $ astarithmetic image.fits -h1        set-in \
                          5 4 in filter-mean    set-filtered \
                          filtered in isblank nan where \
                          --output=out.fits

‘filter-median’
     Apply median filtering
     (https://en.wikipedia.org/wiki/Median_filter) on the input dataset.
     This is very similar to ‘filter-mean’, except that instead of the
     mean value of the box pixels, the median value is used to replace a
     pixel's value.  For more on how to use this operator, please see
     ‘filter-mean’.

     The median is less susceptible to outliers compared to the mean.
     As a result, after median filtering, the pixel values will be more
     discontinuous than mean filtering.

‘filter-minimum’
‘filter-maximum’
     Apply a minimum/maximum filter to the input.  This is very similar
     to ‘filter-mean’, except that instead of the mean value of the box
     pixels, the minimum/maximum value is used to replace a pixel's
     value.  For more on how to use this operator, please see
     ‘filter-mean’.

‘filter-sigclip-mean’
     Apply a $\sigma$-clipped mean filtering onto the input dataset.
     This is very similar to ‘filter-mean’, except that all outliers
     (identified by the $\sigma$-clipping algorithm) have been removed,
     see *note Sigma clipping:: for more on the basics of this
     algorithm.  As described there, two extra input parameters are
     necessary for $\sigma$-clipping: the multiple of $\sigma$ and the
     termination criteria.  ‘filter-sigclip-mean’ therefore needs to pop
     two other operands from the stack after the dimensions of the box.

     For example, the line below uses the same box size as the example
     of ‘filter-mean’.  However, all elements in the box that are
     iteratively beyond $3\sigma$ of the distribution's median are
     removed from the final calculation of the mean until the change in
     $\sigma$ is less than $0.2$.

          $ astarithmetic 3 0.2 5 4 image.fits filter-sigclip-mean

     The median (which needs a sorted dataset) is necessary for
     $\sigma$-clipping, therefore ‘filter-sigclip-mean’ can be
     significantly slower than ‘filter-mean’.  However, if there are
     strong outliers in the dataset that you want to ignore (for
     example, emission lines on a spectrum when finding the continuum),
     this is a much better solution.  Note that $\sigma$-clipping is
     easily biased by outliers, so ‘filter-madclip-mean’ is more robust.

‘filter-sigclip-std’
‘filter-sigclip-mad’
‘filter-sigclip-median’
     Apply a $\sigma$-clipped STD/MAD/median filtering onto the input
     dataset.  This operator and its necessary operands are almost
     identical to ‘filter-sigclip-mean’, except that after
     $\sigma$-clipping, the median value (which is less affected by
     outliers than the mean) is added back to the stack.  Note that
     $\sigma$-clipping is easily biased by outliers, so
     ‘filter-madclip-median’ is more robust.

‘filter-madclip-mean’
     Apply a MAD-clipped mean filtering onto the input dataset.  MAD
     stands for Median Absolute Deviation, this is much more robust to
     outliers than the $\sigma$-clipping.  This operator is called in a
     similar way to ‘filter-madclip-mean’, for example in the command
     below, we use $4\times$MAD clipping with a termination criteria of
     0.01 on an image:
          $ astarithmetic 4 0.01 5 5 image.fits filter-madclip-mean

‘filter-madclip-std’
‘filter-madclip-mad’
‘filter-madclip-median’
     Apply a MAD-clipped STD/MAD/median filtering onto the input
     dataset; see the description of ‘filter-madclip-mean’ for more.


File: gnuastro.info,  Node: Pooling operators,  Next: Interpolation operators,  Prev: Filtering operators,  Up: Arithmetic operators

6.2.4.9 Pooling operators
.........................

Pooling is one way of reducing the complexity of the input image by
grouping multiple input pixels into one output pixel (using any
statistical measure).  As a result, the output image has fewer pixels
(less complexity).  In Computer Vision, Pooling is commonly used in
Convolutional Neural Networks
(https://en.wikipedia.org/wiki/Convolutional_neural_network) (CNNs).

   In pooling, the inputs are an image (e.g., a FITS file) and a square
window pixel size that is known as a pooling window.  The window has to
be smaller than the input's number of pixels in both dimensions and its
width is called the "pool size".  The pooling window starts at the
top-left corner pixel of the input and calculates statistical operations
on the pixels that overlap with it.  It slides forward by the "stride"
pixels, moving over all pixels in the input from the top-left corner to
the bottom-right corner, and repeats the same calculation for the
overlapping pixels in each position.

   Usually, the stride (or spacing between the windows as they slide
over the input) is equal to the window-size.  In other words, in
pooling, the separate "windows" do not overlap with each other on the
input.  However, you can choose any size for the stride.  Remember this,
It's crucial to ensure that the stride size is less than the pool size.
If not, some pixels may be missed during the pooling process.  Therefore
there are two major differences with *note Spatial domain convolution::
or *note Filtering operators::, but pooling has some similarities to the
*note Warp::.
   • In convolution or filtering the input and output sizes are the
     same.  However, when the stride is larger than 1 then, the output
     of pooling must have fewer pixels.
   • In convolution or filters, the kernels slide over the input in a
     pixel-by-pixel manner.  As a result, the same pixel's value will be
     used in many of the output pixels.  However, in pooling each input
     pixel may be only used in a single output pixel (if the stride and
     the pool size are the same).
   • Special cases of Warping an image are similar to pooling.  For
     example calling ‘pool-sum’ with pool size of 2 will give the same
     pixel values (except the outer edges) as giving the same image to
     ‘astwarp’ with ‘--scale=1/2 --centeroncorner’.  However, warping
     will only provide the sum of the input pixels, there is no easy way
     to generically define something like ‘pool-max’ in Warp (which is
     far more general than pooling).  Also, due to its generic features
     (for example for non-linear warps), Warp is slower than the
     ‘pool-max’ that is introduced here.

*No WCS in output:* As of Gnuastro 0.23.119-586f, the output of pooling
will not contain WCS information (primarily due to a lack of time by
developers).  Please inform us of your interest in having it, by
contacting us at ‘bug-gnuastro@gnu.org’.  If you need ‘pool-sum’, you
can use *note Warp:: (which also modifies the WCS, see note above).

   If the width or height of input is not divisible by the stride size,
the pool window will go beyond the input pixel grid.  In this case, the
window pixels that do not overlap with the input are given a blank value
(and thus ignored in the calculation of the desired statistical
operation).

   The simple ASCII figure below shows the pooling operation where the
input is a $3\times3$ pixel image with a pool size of 2 pixels.  In the
center of the second row, you see the intermediate input matrix that
highlights how the input and output pixels relate with each other.
Since the input is $3\times3$ and we have a stride size of 2, as
mentioned above blank pseudo-pixels are added with a value of ‘B’ (for
blank).

             Pool window:                             Input:
             +-----------+                           +-------------+
             |     |     |                           | 10   12   9 |
             | _ _ | _ _ |___________________________| 31   4    1 |
             |     |     |       ||          ||      | 16   5    8 |
             |     |     |       ||          ||      +-------------+
             +-----------+       ||          ||
         The pooling window 2*2  ||          ||
                stride 2         \/          \/
                             +---------------------+
                             |/ 10   12\|/ 9    B \|
                             |          |          |
       +-------+  pool-min   |\ 31   4 /|\ 1    B /|   pool-max  +-------+
       | 4   1 |   /------   |---------------------|   ------\   |31   9 |
       | 5   8 |   \------   |/ 16   5 \|/ 8    B \|   ------/   |16   8 |
       +-------+             |          |          |             +-------+
                             |\ B    B /.\ B    B /|
                             +---------------------+

   The choice of the statistic to use depends on the specific use case,
the characteristics of the input data, and the desired output.  Each
statistic has its advantages and disadvantages and the choice of which
to use should be informed by the specific needs of the problem at hand.
Below, the various pool operators of arithmetic are listed:

‘pool-max’
     Apply max-pooling on the input dataset.  This operator takes three
     operands: the first popped operand is the stride and the second is
     the width of the square pooling window (which should be a single
     integer).  Also, The third operand should be the input image.
     Within the pooling window, this operator will place the largest
     value in the output pixel (any blank pixels will be ignored).

     See the ASCII diagram above for a demonstration of how max-pooling
     works.  Here is an example of using this operator:

          $ astarithmetic image.fits 2 2 pool-max

     Max-pooling retains the largest value of the input window in the
     output, so the returned image is sharper where you have strong
     signal-to-noise ratio and more noisy in regions with no significant
     signal (only noise).  It is therefore useful when the background of
     the image is dark and we are interested in only the highest
     signal-to-noise ratio regions of the image.

‘pool-min’
     Apply min-pooling on the input dataset.  This operator takes three
     operands: the first popped operand is the stride and the second is
     the width of the square pooling window (which should be a single
     integer).  Also, The third operand should be the input image.
     Except the used statistical measurement, this operator is similar
     to ‘pool-max’, see the description there for more.

     Min-pooling is mostly used when the image has a high
     signal-to-noise ratio and a light background: min-pooling will
     select darker (lower-valued) pixels.  For low signal-to-noise
     regions, this operator will increase the noise level (similar to
     the maximum, the scatter in the minimum is very strong).

‘pool-sum’
     Apply sum-pooling to the input dataset.  This operator takes three
     operands: the first popped operand is the stride and the second is
     the width of the square pooling window (which should be a single
     integer).  Also, The third operand should be the input image.
     Except the used statistical measurement, this operator is similar
     to ‘pool-max’, see the description there for more.

     Sum-pooling will increase the signal-to-noise ratio at the cost of
     having a smoother output (less resolution).

‘pool-mean’
     Apply mean pooling on the input dataset.  This operator takes three
     operands: the first popped operand is the stride and the second is
     the width of the square pooling window (which should be a single
     integer).  Also, The third operand should be the input image.
     Except the used statistical measurement, this operator is similar
     to ‘pool-max’, see the description there for more.

     The mean pooling method smooths out the image and hence the sharp
     features may not be identified when this pooling method is used.
     This therefore preserves more information than max-pooling, but may
     also reduces the effect of the most prominent pixels.  Mean is
     often used where a more accurate representation of the input is
     required.

‘pool-median’
     Apply median pooling on the input dataset.  This operator takes
     three operands: the first popped operand is the stride and the
     second is the width of the square pooling window (which should be a
     single integer).  Also, The third operand should be the input
     image.  Except the used statistical measurement, this operator is
     similar to ‘pool-max’, see the description there for more.

     In general, the mean is mathematically easier to interpret and more
     susceptible to outliers, while the median outputs as being less
     subject to the influence of outliers compared to the mean so we
     have a smoother image.  This is therefore better for low
     signal-to-ratio (noisy) features and extended features (where you
     don't want a single high or low valued pixel to affect the output).


File: gnuastro.info,  Node: Interpolation operators,  Next: Dimensionality changing operators,  Prev: Pooling operators,  Up: Arithmetic operators

6.2.4.10 Interpolation operators
................................

Interpolation is the process of removing blank pixels from a dataset (by
giving them a value based on the non-blank neighbors).

‘interpolate-medianngb’
     Interpolate the blank elements of the second popped operand with
     the median of nearest non-blank neighbors to each.  The number of
     the nearest non-blank neighbors used to calculate the median is
     given by the first popped operand.

     The distance of the nearest non-blank neighbors is irrelevant in
     this interpolation.  The neighbors of each blank pixel will be
     parsed in expanding circular rings (for 2D images) or spherical
     surfaces (for 3D cube) and each non-blank element over them is
     stored in memory.  When the requested number of non-blank neighbors
     have been found, their median is used to replace that blank
     element.  For example, the line below replaces each blank element
     with the median of the nearest 5 pixels.

          $ astarithmetic image.fits 5 interpolate-medianngb

     When you want to interpolate blank regions and you want each blank
     region to have a fixed value (for example, the centers of saturated
     stars) this operator is not good.  Because the pixels used to
     interpolate various parts of the region differ.  For such
     scenarios, you may use ‘interpolate-maxofregion’ or
     ‘interpolate-inofregion’ (described below).

‘interpolate-meanngb’
     Similar to ‘interpolate-medianngb’, but will fill the blank values
     of the dataset with the mean value of the requested number of
     nearest neighbors.

‘interpolate-minngb’
     Similar to ‘interpolate-medianngb’, but will fill the blank values
     of the dataset with the minimum value of the requested number of
     nearest neighbors.

‘interpolate-maxngb’
     Similar to ‘interpolate-medianngb’, but will fill the blank values
     of the dataset with the maximum value of the requested number of
     nearest neighbors.  One useful implementation of this operator is
     to fill the saturated pixels of stars in images.

‘interpolate-minofregion’

     Interpolate all blank regions (consisting of many blank pixels that
     are touching) of the third popped operand with the minimum value of
     the pixels that are touching that region (a single value).  The
     second popped operand is the inner-connectivity (defining the
     actual region: which touching pixels are considered as part of the
     region to give a single value).  The first popped operand is the
     outer-connectivity which defines the pixels outside the blank
     region which are used.  For more on the definition of connectivity,
     see ‘connected-components’).

     For example, with the command below all the 2-connected blank
     regions of ‘image.fits’ will be filled with the minimum of the
     1-connected, non-blank neighbors of the region.  The 'image.fits'
     is an image (2D dataset), so a 2 connectivity means that the
     independent blank regions are defined by 8-connected neighbors: two
     pixels touching on the corner (a point) are considered as one
     region.  If inner-connectivity was 1, the regions would be defined
     by 4-connectivity: blank regions would be defined only by pixels
     touch touch along one full edge length.  Conversely, for the
     non-blank pixels "touching" the outer parts of the blank regions
     (to find the minimum), the command below assumes 1-connectivity.

          $ astarithmetic image.fits 2 1 interpolate-minofregion

‘interpolate-maxofregion’
     Similar to ‘interpolate-minofregion’, but the maximum is used to
     fill the blank regions.

     This operator can be useful in filling saturated pixels in stars
     for example.  Recall that the ‘interpolate-maxngb’ operator looks
     for the maximum value with a given number of neighboring pixels and
     is more useful in small noisy regions.  Therefore as the blank
     regions become larger, ‘interpolate-maxngb’ can cause a
     fragmentation in the connected blank region because the nearest
     neighbor to one part of the blank region, may not fall within the
     pixels searched for the other regions.  With this option, the size
     of the blank region is irrelevant: all the pixels bordering the
     blank region are parsed and their maximum value is used for the
     whole region.


File: gnuastro.info,  Node: Dimensionality changing operators,  Next: Conditional operators,  Prev: Interpolation operators,  Up: Arithmetic operators

6.2.4.11 Dimensionality changing operators
..........................................

Through these operators you can change the dimensions of the output
through certain statistics on the dimensions that should be removed.
For example, let's assume you have a 3D data cube that has 300 by 300
pixels in the RA and Dec dimensions (first two dimensions), and 3600
slices along the wavelength (third dimension), so the whole cube is
$300\times300\times3600$ voxels (volume elements).  To create a
narrow-band image that only contains 100 slices around a certain
wavelength, you can crop that section (using *note Crop::), giving you a
$300\times300\times100$ cube.  You can now use the ‘collapse-sum’
operator below to "collapse" all the 100 slices into one 2D image that
has $300\times300$ pixels.  Every pixel in this 2D image will have the
flux of the sum of the 100 slices.

‘to-1d’
     Convert the input operand into a 1D array; irrespective of the
     number of dimensions it has.  This operator only takes a single
     operand (the input array) and just updates the metadata.  Therefore
     it does not change the layout of the array contents in memory and
     is very fast.

     If no further operation is requested on the 1D array, recall that
     Arithmetic will write a 1D array as a table column by default.  In
     case you want the output to be saved as a 1D image, or to see it on
     the standard output, please use the ‘--onedasimage’ or
     ‘--onedonstdout’ options respectively (see *note Invoking
     astarithmetic::).

     This operator is useful in scenarios where after some operations on
     a 2D image or 3D cube, the dimensionality is no longer relevant for
     you and you just care about the values.  In the example below, we
     will first make a simple 2D image from a plain-text file, then
     convert it to a 1D array:

          ## Contents of 'a.txt' to start with.
          $ cat a.txt
          # Image 1: DEMO [counts, uint8] An example image
          1 2 3
          4 5 6
          7 8 9

          ## Convert the text image into a FITS image.
          $ astconvertt a.txt -o a.fits

          ## Convert it into a table column (1D):
          $ astarithmetic a.fits to-1d -o table.fits

          ## Convert it into a 1D image:
          $ astarithmetic a.fits to-1d -o table.fits --onedasimage

     A more real-world example would be the following: assume you want
     to "flatten" two images into a single 1D array (as commonly done in
     convolutional neural networks, or CNNs(1)).  First, we show the
     contents of a new $2\times2$ image in plain-text image, then
     convert it to a 2D FITS image (‘b.fits’).  We will then use
     arithmetic to make both ‘a.fits’ (from the example above) and
     ‘b.fits’ into a 1D array and stitch them together into a single 1D
     image with one call to Arithmetic.  For a description of the
     ‘stitch’ operator, see below (same section).

          ## Contents of 'b.txt':
          $ cat b.txt
          # Image 1: DEMO [counts, uint8] An example image
          10 11
          12 13

          ## Convert the text image into a FITS image.
          $ astconvertt b.txt -o b.fits

          # Flatten the two images into a single 1D image:
          $ astarithmetic a.fits to-1d b.fits to-1d 2 1 stitch -g1 \
                          --onedonstdout --quiet
          1
          2
          3
          4
          5
          6
          7
          8
          9
          10
          11
          12
          13

‘stitch’
     Stitch (connect) any number of given images together along the
     given dimension.  The output has the same number of dimensions as
     the input, but the number of pixels along the requested dimension
     will be different from the inputs.  The ‘stitch’ operator takes at
     least three operands:
        • The first popped operand (placed just before ‘stitch’) is the
          direction (dimension) that the images should be stitched
          along.  The first FITS dimension is along the horizontal,
          therefore a value of ‘1’ will stitch them horizontally.
          Similarly, giving a value of ‘2’ will result in a vertical
          stitch.

        • The second popped operand is the number of images that should
          be stitched.

        • Depending on the value given to the second popped operand,
          ‘stitch’ will pop the given number of datasets from the stack
          and stitch them along the given dimension.  The popped images
          have to have the same number of pixels along the other
          dimension.  The order of the stitching is defined by how they
          are placed in the command-line, not how they are popped (after
          being popped, they are placed in a list in the same order).

     For example, in the commands below, we will first crop out fixed
     sized regions of $100\times300$ pixels of a larger image
     (‘large.fits’) first.  In the first call of Arithmetic below, we
     will stitch the bottom set of crops together along the first
     (horizontal) axis.  In the second Arithmetic call, we will stitch
     all 6 along both dimensions.

          ## Crop the fixed-size regions of a larger image ('-O' is the
          ## short form of the '--mode' option).
          $ astcrop large.fits -Oimg --section=1:100,1:300     -oa.fits
          $ astcrop large.fits -Oimg --section=101:200,1:300   -ob.fits
          $ astcrop large.fits -Oimg --section=201:300,1:300   -oc.fits
          $ astcrop large.fits -Oimg --section=1:100,301:600   -od.fits
          $ astcrop large.fits -Oimg --section=101:200,301:600 -oe.fits
          $ astcrop large.fits -Oimg --section=201:300,301:600 -of.fits

          ## Stitch the bottom three crops into one image.
          $ astarithmetic a.fits b.fits c.fits 3 1 stitch -obottom.fits

          # Stitch all the 6 crops along both dimensions
          $ astarithmetic a.fits b.fits c.fits 3 1 stitch \
                          d.fits e.fits f.fits 3 1 stitch \
                          2 2 stitch -g1 -oall.fits

     The start of the last command is like the one before it (stitching
     the bottom three crops along the first FITS dimension, producing a
     $300\times300$ image).  Later in the same command, we then stitch
     the top three crops horizontally (again, into a $300\times300$
     image) This leaves the the two $300\times300$ images on the stack
     (see *note Reverse polish notation::).  We finally stitch those two
     along the second (vertical) dimension.  This operator is therefore
     useful in scenarios like placing the CCD amplifiers into one image.

‘trim’
     Trim all blank elements from the outer edges of the input operand
     (it only takes a single operand).  For example see the commands
     below using Table's *note Column arithmetic:::

          $  cat table.txt
          nan
          nan
          nan
          3
          4
          nan
          5
          6
          nan

          $ asttable table.txt -Y -c'arith $1 trim'
          3.000000
          4.000000
          nan
          5.000000
          6.000000

     Similarly, on 2D images or 3D cubes, all outer rows/columns or
     slices that are fully blank get "trim"ed with this operator.  This
     is therefore a very useful operator for extracting a certain
     feature within your dataset.

     For example, let's assume that you have set *note NoiseChisel:: and
     *note Segment:: on an image to extract all clumps and objects.
     With the command below on Segment's output, you will have a smaller
     image that only contains the sky-subtracted input pixels
     corresponding to object 263.

          $ astarithmetic seg.fits -hINPUT-NO-SKY seg.fits -hOBJECTS \
                          263 ne nan where trim --output=obj-263.fits

‘add-dimension-slow’
     Build a higher-dimensional dataset from all the input datasets
     stacked after one another (along the slowest dimension).  The first
     popped operand has to be a single number.  It is used by the
     operator to know how many operands it should pop from the stack
     (and the size of the output in the new dimension).  The rest of the
     operands must have the same size and numerical data type.  This
     operator currently only works for 2D input operands, please contact
     us if you want inputs to have different dimensions.

     The output's WCS (which should have a different dimensionality
     compared to the inputs) can be read from another file with the
     ‘--wcsfile’ option.  If no file is specified for the WCS, the first
     dataset's WCS will be used, you can later add/change the necessary
     WCS keywords with the FITS keyword modification features of the
     Fits program (see *note Fits::).

     If your datasets do not have the same type, you can use the type
     transformation operators of Arithmetic that are discussed below.
     Just beware of overflow if you are transforming to a smaller type,
     see *note Numeric data types::.

     For example, let's assume you have 3 two-dimensional images
     ‘a.fits’, ‘b.fits’ and ‘c.fits’ (each with $200\times100$ pixels).
     You can construct a 3D data cube with $200\times100\times3$ voxels
     (volume-pixels) using the command below:

          $ astarithmetic a.fits b.fits c.fits 3 add-dimension-slow

‘add-dimension-fast’
     Similar to ‘add-dimension-slow’ but along the fastest dimension.
     This operator currently only works for 1D input operands, please
     contact us if you want inputs to have different dimensions.

     For example, let's assume you have 3 one-dimensional datasets, each
     with 100 elements.  With this operator, you can construct a
     $3\times100$ pixel FITS image that has 3 pixels along the
     horizontal and 5 pixels along the vertical.

‘collapse-sum’
     Collapse the given dataset (second popped operand), by summing all
     elements along the first popped operand (a dimension in FITS
     standard: counting from one, from fastest dimension).  The returned
     dataset has one dimension less compared to the input.

     The output will have a double-precision floating point type
     irrespective of the input dataset's type.  Doing the operation in
     double-precision (64-bit) floating point will help the collapse
     (summation) be affected less by floating point errors.  But
     afterwards, single-precision floating points are usually enough in
     real (noisy) datasets.  So depending on the type of the input and
     its nature, it is recommended to use one of the type conversion
     operators on the returned dataset.

     If any WCS is present, the returned dataset will also lack the
     respective dimension in its WCS matrix.  Therefore, when the WCS is
     important for later processing, be sure that the input is aligned
     with the respective axes: all non-diagonal elements in the WCS
     matrix are zero.

     One common application of this operator is the creation of
     synthetic broad-band or narrow-band 2D images from 3D data cubes.
     For example, integral field unit (IFU) data products that have two
     spatial dimensions (first two FITS dimensions) and one spectral
     dimension (third FITS dimension).  The command below will collapse
     the whole third dimension into a 2D array the size of the first two
     dimensions, and then convert the output to single-precision
     floating point (as discussed above).

          $ astarithmetic cube.fits 3 collapse-sum float32

‘collapse-min’
‘collapse-max’
‘collapse-mean’
‘collapse-median’
     Similar to ‘collapse-sum’, but the returned dataset will be the
     desired statistic along the collapsed dimension, not the sum.

‘collapse-madclip-fill-mad’
‘collapse-madclip-fill-std’
‘collapse-madclip-fill-mean’
‘collapse-madclip-fill-median’
‘collapse-madclip-fill-number’
     Collapse the input dataset (fourth popped operand) along the FITS
     dimension given as the first popped operand by calculating the
     desired statistic after median absolute deviation (MAD) filled
     re-clipping.  The MAD-clipping parameters (namely, the multiple of
     sigma and termination criteria) are read as the third and second
     popped operands respectively.

     This is the most robust method to reject outliers; for more on
     filled re-clipping and its advantages, see *note Contiguous
     outliers::.  For a more general tutorial on rejecting outliers, see
     *note Clipping outliers::.  If you have not done this tutorial yet,
     we recommend you to take an hour or so and go through that tutorial
     for optimal understanding and results.

     When more than 95$\%$ of the area of an operand is masked, the full
     operand will be masked.  This is necessary in like this: one of
     your inputs has many outliers (for example it is much more noisy
     than the rest or its sky level has not been subtracted properly).
     Because this operator fills holes between outlying pixels, most of
     the area of the input will be masked, but the thin edges (where
     there are no "holes") will remain, causing different statistics in
     those thin edges of that input in your final coadd.  Through this
     mask coverage fraction (which is currently hard-coded(2)), we
     ensure that such thin edges do not cause artifacts in the final
     coadd.

     For example, with the command below, the pixels of the input 2
     dimensional ‘image.fits’ will be collapsed to a single dimension
     output.  The first popped operand is ‘2’, so it will collapse all
     the pixels that are vertically on top of each other.  Such that the
     output will have the same number of pixels as the horizontal axis
     of the input.  During the collapsing, all pixels that are more than
     $3\sigma$ (third popped operand) are rejected, and the clipping
     will continue until the standard deviation changes less than $0.2$
     between clips.  Finally the ‘counter’ operator is used to have a
     two-column table with the first one being a simple counter starting
     from one (see *note Size and position operators::).

          $ astarithmetic image.fits 3 0.2 2 collapse-sigclip-mean \
                          counter --output=collapsed-vertical.fits

     *Printing output of collapse in plain-text:* the default datatype
     of ‘collapse-sigclip-mean’ is 32-bit floating point.  This is
     sufficient for any observed astronomical data.  However, if you
     request a plain-text output, or decide to print/view the output as
     plain-text on the standard output, the full set of decimals may not
     be printed in some situations.  This can lead to apparently
     discrete values in the output of this operator when viewed in
     plain-text!  The FITS format is always superior (since it stores
     the value in binary, therefore not having the problem above).  But
     if you are forced to save the output in plain-text, use the
     ‘float64’ operator after this to change the type to 64-bit floating
     point (which will print more decimals).

‘collapse-madclip-mad’
‘collapse-madclip-std’
‘collapse-madclip-mean’
‘collapse-madclip-median’
‘collapse-madclip-number’
     Collapse the input dataset (fourth popped operand) along the FITS
     dimension given as the first popped operand by calculating the
     desired statistic after median absolute deviation (MAD) clipping.
     This operator is called similarly to the ‘collapse-madclip-fill-*’
     operators, see the description there for more.

‘collapse-sigclip-fill-mad’
‘collapse-sigclip-fill-std’
‘collapse-sigclip-fill-mean’
‘collapse-sigclip-fill-median’
‘collapse-sigclip-fill-number’
     Collapse the input dataset (fourth popped operand) along the FITS
     dimension given as the first popped operand by calculating the
     desired statistic after filled $\sigma$ re-clipping.  This operator
     is called similarly to the ‘collapse-madclip-fill-*’ operators, see
     the description there for more.

‘collapse-sigclip-mad’
‘collapse-sigclip-std’
‘collapse-sigclip-mean’
‘collapse-sigclip-median’
‘collapse-sigclip-number’
     Collapse the input dataset (fourth popped operand) along the FITS
     dimension given as the first popped operand by calculating the
     desired statistic after $\sigma$-clipping.  This operator is called
     similarly to the ‘collapse-madclip-fill-*’ operators, see the
     description there for more.

   ---------- Footnotes ----------

   (1) <https://en.wikipedia.org/wiki/Convolutional_neural_network>

   (2) Please get in touch with us at ‘bug-gnuastro@gnu.org’ if you
notice this problem and feel the fraction needs to be lowered (or
generally to be set in each run).


File: gnuastro.info,  Node: Conditional operators,  Next: Mathematical morphology operators,  Prev: Dimensionality changing operators,  Up: Arithmetic operators

6.2.4.12 Conditional operators
..............................

Conditional operators take two inputs and return a binary output that
can only have two values 0 (for pixels where the condition was false) or
1 (for the pixels where the condition was true).  Because of the binary
(2-valued) nature of their outputs, the output is therefore stored in an
‘unsigned char’ data type (see *note Numeric data types::) to speed up
process and take less space in your storage.  There are two exceptions
to the general features above: ‘isblank’ only takes one input, and
‘where’ takes three, while not returning a binary output, see their
description for more.

‘lt’
     Less than: creates a binary output (values either 0 or 1) where
     each pixel will be 1 if the second popped operand is smaller than
     the first popped operand and 0 otherwise.  If both operands are
     images, then all the pixels will be compared with their
     counterparts in the other image.

     For example, the pixels in the output of the command below will
     have a value of 1 (true) if their value in ‘image1.fits’ is less
     than their value in ‘image2.fits’.  Otherwise, their value will be
     0 (false).
          $ astarithmetic image1.fits image2.fits lt
     If only one operand is an image, then all the pixels will be
     compared with the single value (number) of the other operand.  For
     example:
          $ astarithmetic image1.fits 1000 lt
     Finally if both are numbers, then the output is also just one
     number (0 or 1).
          $ astarithmetic 4 5 lt

‘le’
     Less or equal: similar to ‘lt’ ('less than' operator), but
     returning 1 when the second popped operand is smaller or equal to
     the first.  For example
          $ astarithmetic image1.fits 1000 le

‘gt’
     Greater than: similar to ‘lt’ ('less than' operator), but returning
     1 when the second popped operand is greater than the first.  For
     example
          $ astarithmetic image1.fits 1000 gt

‘ge’
     Greater or equal: similar to ‘lt’ ('less than' operator), but
     returning 1 when the second popped operand is larger or equal to
     the first.  For example
          $ astarithmetic image1.fits 1000 ge

‘eq’
     Equality: similar to ‘lt’ ('less than' operator), but returning 1
     when the two popped operands are equal (to double precision
     floating point accuracy).
          $ astarithmetic image1.fits 1000 eq

‘ne’
     Non-Equality: similar to ‘lt’ ('less than' operator), but returning
     1 when the two popped operands are _not_ equal (to double precision
     floating point accuracy).
          $ astarithmetic image1.fits 1000 ne

‘and’
     Logical AND: returns 1 if both operands have a non-zero value, and
     returns 0 if either operand is zero.  Both operands have to be the
     same kind: either both images or both numbers and it mostly makes
     meaningful values when the inputs are binary (with pixel values of
     0 or 1).
          $ astarithmetic image1.fits image2.fits -g1 and

     For example, if you only want to see which pixels in an image have
     a value _between_ 50 (greater equal, or inclusive) and 200 (less
     than, or exclusive), you can use this command:
          $ astarithmetic image.fits set-i i 50 ge i 200 lt and

‘or’
     Logical OR: returns 1 if either one of the operands is non-zero and
     0 only when both operators are zero.  Both operands have to be the
     same kind: either both images or both numbers.  The usage is
     similar to ‘and’.

     For example, if you only want to see which pixels in an image have
     a value _outside of_ -100 (greater equal, or inclusive) and 200
     (less than, or exclusive), you can use this command:
          $ astarithmetic image.fits set-i i -100 lt i 200 ge or

‘not’
     Logical NOT: returns 1 when the operand is 0 and 0 when the operand
     is non-zero.  The operand can be an image or number, for an image,
     it is applied to each pixel separately.  For example, if you want
     to know which pixels are not blank (and assuming that we didn't
     have the ‘isnotblank’ operator), you can use this ‘not’ operator on
     the output of the ‘isblank’ operator described below:
          $ astarithmetic image.fits isblank not

‘isblank’
     Test each pixel for being a blank value (see *note Blank pixels::).
     This is a conditional operator: the output has the same size and
     dimensions as the input, but has an unsigned 8-bit integer type
     with two possible values: either 1 (for a pixel that was blank) or
     0 (for a pixel that was not blank).  See the description of ‘lt’
     operator above).  The difference is that it only needs one operand.
     For example:
          $ astarithmetic image.fits isblank
     Because of the definition of a blank pixel, a blank value is not
     even equal to itself, so you cannot use the equal operator above to
     select blank pixels.  See the "Blank pixels" box below for more on
     Blank pixels in Arithmetic.

     In case you want to set non-blank pixels to an output pixel value
     of 1, it is better to use ‘isnotblank’ instead of '‘isblank not’'
     (for more, see the description of ‘isnotblank’).

‘isnotblank’
     The inverse of the ‘isblank’ operator above (see that description
     for more).  Therefore, if a pixel has a blank value, the output of
     this operator will have a 0 value for it.  This operator is
     therefore similar to running '‘isblank not’', but slightly more
     efficient (won't need the intermediate product of two operators).

‘where’
     Change the input (pixel) value _where_/if a certain condition
     holds.  The conditional operators above can be used to define the
     condition.  Three operands are required for ‘where’.  The input
     format is demonstrated in this simplified example:

          $ astarithmetic modify.fits binary.fits if-true.fits where

     The value of any pixel in ‘modify.fits’ that corresponds to a
     non-zero _and_ non-blank pixel of ‘binary.fits’ will be changed to
     the value of the same pixel in ‘if-true.fits’ (this may also be a
     number).  The 3rd and 2nd popped operands (‘modify.fits’ and
     ‘binary.fits’ respectively, see *note Reverse polish notation::)
     have to have the same dimensions/size.  ‘if-true.fits’ can be
     either a number, or have the same dimension/size as the other two.

     The 2nd popped operand (‘binary.fits’) has to have ‘uint8’ (or
     ‘unsigned char’ in standard C) type (see *note Numeric data
     types::).  It is treated as a binary dataset (with only two values:
     zero and non-zero, hence the name ‘binary.fits’ in this example).
     However, commonly you will not be dealing with an actual FITS file
     of a condition/binary image.  You will probably define the
     condition in the same run based on some other reference image and
     use the conditional and logical operators above to make a
     true/false (or one/zero) image for you internally.  For example,
     the case below:

          $ astarithmetic in.fits reference.fits 100 gt new.fits where

     In the example above, any of the ‘in.fits’ pixels that has a value
     in ‘reference.fits’ greater than ‘100’, will be replaced with the
     corresponding pixel in ‘new.fits’.  Effectively the ‘reference.fits
     100 gt’ part created the condition/binary image which was added to
     the stack (in memory) and later used by ‘where’.  The command above
     is thus equivalent to these two commands:

          $ astarithmetic reference.fits 100 gt --output=binary.fits
          $ astarithmetic in.fits binary.fits new.fits where

     Finally, the input operands are read and used independently, so you
     can use the same file more than once as any of the operands.

     When the 1st popped operand to ‘where’ (‘if-true.fits’) is a single
     number, it may be a NaN value (or any blank value, depending on its
     type) like the example below (see *note Blank pixels::).  When the
     number is blank, it will be converted to the blank value of the
     type of the 3rd popped operand (‘in.fits’).  Hence, in the example
     below, all the pixels in ‘reference.fits’ that have a value greater
     than 100, will become blank in the natural data type of ‘in.fits’
     (even though NaN values are only defined for floating point types).

          $ astarithmetic in.fits reference.fits 100 gt nan where


File: gnuastro.info,  Node: Mathematical morphology operators,  Next: Bitwise operators,  Prev: Conditional operators,  Up: Arithmetic operators

6.2.4.13 Mathematical morphology operators
..........................................

From Wikipedia: "Mathematical morphology (MM) is a theory and technique
for the analysis and processing of geometrical structures, based on set
theory, lattice theory, topology, and random functions.  MM is most
commonly applied to digital images".  In theory it extends a very large
body of research and methods in image processing, but currently in
Gnuastro it mainly applies to images that are binary (only have a value
of 0 or 1).  For example, you have applied the greater-than operator
(‘gt’, see *note Conditional operators::) to select all pixels in your
image that are larger than a value of 100.  But they will all have a
value of 1, and you want to separate the various groups of pixels that
are connected (for example, peaks of stars in your image).  With the
‘connected-components’ operator, you can give each connected region of
the output of ‘gt’ a separate integer label.

‘erode’
     Erode the foreground pixels (with value ‘1’) of the input dataset
     (second popped operand).  The first popped operand is the
     connectivity (see description in ‘connected-components’).  Erosion
     is simply a flipping of all foreground pixels (with value ‘1’) to
     background (with value ‘0’) that are "touching" background pixels.
     "Touching" is defined by the connectivity.

     In effect, this operator "carves off" the outer borders of the
     foreground, making them thinner.  This operator assumes a binary
     dataset (all pixels are ‘0’ or ‘1’).  For example, imagine that you
     have an astronomical image with a mean/sky value of 0 units and a
     standard deviation ($\sigma$) of 100 units and many galaxies in it.
     With the first command below, you can apply a threshold of
     $2\sigma$ on the image (by only keeping pixels that are greater
     than 200 using the ‘gt’ operator).  The output of thresholding the
     image is a binary image (each pixel is either smaller or equal to
     the threshold or larger than it).  You can then erode the binary
     image with the second command below to remove very small false
     positives (one or two pixel peaks).
          $ astarithmetic image.fits 100 gt -obinary.fits
          $ astarithmetic binary.fits 2 erode -oout.fits

     In fact, you can merge these operations into one command thanks to
     the reverse polish notation (see *note Reverse polish notation::):
          $ astarithmetic image.fits 100 gt 2 erode -oout.fits

     To see the effect of connectivity, try this:
          $ astarithmetic image.fits 100 gt 1 erode -oout-con-1.fits

‘dilate’
     Dilate the foreground pixels (with value ‘1’) of the binary input
     dataset (second popped operand).  The first popped operand is the
     connectivity (see description in ‘connected-components’).  Dilation
     is simply a flipping of all background pixels (with value ‘0’) to
     foreground (with value ‘1’) that are "touching" foreground pixels.
     "Touching" is defined by the connectivity.  In effect, this expands
     the outer borders of the foreground.  This operator assumes a
     binary dataset (all pixels are ‘0’ and ‘1’).  The usage is similar
     to ‘erode’, for example:
          $ astarithmetic binary.fits 2 dilate -oout.fits

‘number-neighbors’
     Return a dataset of the same size as the second popped operand, but
     where each non-zero and non-blank input pixel is replaced with the
     number of its non-zero and non-blank neighbors.  The first popped
     operand is the connectivity (see above) and must be a single-value
     of an integer type.  The dataset is assumed to be binary (having an
     unsigned, 8-bit dataset).

     For example with the command below, you can select all pixels above
     a value of 100 in your image with the "greater-than" or ‘gt’
     operator (see *note Conditional operators::).  Recall that the
     output of all conditional operators is a binary output (having a
     value of 0 or 1).  In the same command, we will then find how many
     neighboring pixels of each pixel (that was originally above the
     threshold) are also above the threshold.

          $ astarithmetic image.fits 100 gt 2 number-neighbors

‘connected-components’
     Find the connected components in the input dataset (second popped
     operand).  The first popped is the connectivity used in the
     connected components algorithm.  The second popped operand is the
     dataset where connected components are to be found.  It is assumed
     to be a binary image (with values of 0 or 1).  It must have an
     8-bit unsigned integer type which is the format produced by
     conditional operators.  This operator will return a labeled dataset
     where the non-zero pixels in the input will be labeled with a
     counter (starting from 1).

     The connectivity is a number between 1 and the number of dimensions
     in the dataset (inclusive).  1 corresponds to the weakest
     (symmetric) connectivity between elements and the number of
     dimensions the strongest.  For example, on a 2D image, a
     connectivity of 1 corresponds to 4-connected neighbors and 2
     corresponds to 8-connected neighbors.

     One example usage of this operator can be the identification of
     regions above a certain threshold, as in the command below.  With
     this command, Arithmetic will first separate all pixels greater
     than 100 into a binary image (where pixels with a value of 1 are
     above that value).  Afterwards, it will label all those that are
     connected.

          $ astarithmetic in.fits 100 gt 2 connected-components

     If your input dataset does not have a binary type, but you know all
     its values are 0 or 1, you can use the ‘uint8’ operator (below) to
     convert it to binary.

‘fill-holes’
     Flip background (0) pixels surrounded by foreground (1) in a binary
     dataset.  This operator takes two operands (similar to
     ‘connected-components’): the second is the binary (0 or 1 valued)
     dataset to fill holes in and the first popped operand is the
     connectivity (to define a hole).  Imagine that in your dataset
     there are some holes with zero value inside the objects with one
     value (for example, the output of the thresholding example of
     ‘erode’) and you want to fill the holes:
          $ astarithmetic binary.fits 2 fill-holes

‘invert’
     Invert an unsigned integer dataset (will not work on other data
     types, see *note Numeric data types::).  This is the only operator
     that ignores blank values (which are set to be the maximum values
     in the unsigned integer types).

     This is useful in cases where the target(s) has(have) been imaged
     in absorption as raw formats (which are unsigned integer types).
     With this option, the maximum value for the given type will be
     subtracted from each pixel value, thus "inverting" the image, so
     the target(s) can be treated as emission.  This can be useful when
     the higher-level analysis methods/tools only work on emission
     (positive skew in the noise, not negative).
          $ astarithmetic image.fits invert


File: gnuastro.info,  Node: Bitwise operators,  Next: Numerical type conversion operators,  Prev: Mathematical morphology operators,  Up: Arithmetic operators

6.2.4.14 Bitwise operators
..........................

Astronomical images are usually stored as an array multi-byte pixels
with different sizes for different precision levels (see *note Numeric
data types::).  For example, images from CCDs are usually in the
unsigned 16-bit integer type (each pixel takes 16 bits, or 2 bytes, of
memory) and fully reduced deep images have a 32-bit floating point type
(each pixel takes 32 bits or 4 bytes).

   On the other hand, during the data reduction, we need to preserve a
lot of meta-data about some pixels.  For example, if a cosmic ray had
hit the pixel during the exposure, or if the pixel was saturated, or is
known to have a problem, or if the optical vignetting is too strong on
it.  A crude solution is to make a new image when checking for each one
of these things and make a binary image where we flag (set to 1) pixels
that satisfy any of these conditions above, and set the rest to zero.
However, processing pipelines sometimes need more than 20 flags to store
important per-pixel meta-data, and recall that the smallest numeric data
type is one byte (or 8 bits, that can store up to 256 different values),
while we only need two values for each flag!  This is a major waste of
storage space!

   A much more optimal solution is to use the bits within each pixel to
store different flags!  In other words, if you have an 8-bit pixel, use
each bit as a flag to mark if a certain condition has happened on a
certain pixel or not.  For example, let's set the following standard
based on the four cases mentioned above: the first bit will show that a
cosmic ray has hit that pixel.  So if a pixel is only affected by cosmic
rays, it will have this sequence of bits (note that the bit-counting
starts from the right): ‘00000001’.  The second bit shows that the pixel
was saturated (‘00000010’), the third bit shows that it has known
problems (‘00000100’) and the fourth bit shows that it was affected by
vignetting (‘00001000’).

   Since each bit is independent, we can thus mark multiple metadata
about that pixel in the actual image, within a single "flag" or "mask"
pixel of a flag or mask image that has the same number of pixels.  For
example, a flag-pixel with the following bits ‘00001001’ shows that it
has been affected by cosmic rays _and_ it has been affected by
vignetting at the same time.  The common data type to store these
flagging pixels are unsigned integer types (see *note Numeric data
types::).  Therefore when you open an unsigned 8-bit flag image in a
viewer like DS9, you will see a single integer in each pixel that
actually has 8 layers of metadata in it!  For example, the integer you
will see for the bit sequences given above will respectively be: $2^0=1$
(for a pixel that only has cosmic ray), $2^1=2$ (for a pixel that was
only saturated), $2^2=4$ (for a pixel that only has known problems),
$2^3=8$ (for a pixel that is only affected by vignetting) and $2^0 + 2^3
= 9$ (for a pixel that has a cosmic ray _and_ was affected by
vignetting).

   You can later use this bit information to mark objects in your final
analysis or to mask certain pixels.  For example, you may want to set
all pixels affected by vignetting to NaN, but can interpolate over
cosmic rays.  You therefore need ways to separate the pixels with a
desired flag(s) from the rest.  It is possible to treat a flag pixel as
a single integer (and try to define certain ranges in value to select
certain flags).  But a much more easier and robust way is to actually
look at each pixel as a sequence of bits (not as a single integer!)  and
use the bitwise operators below for this job.  For more on the theory
behind bitwise operators, see Wikipedia
(https://en.wikipedia.org/wiki/Bitwise_operation).

‘bitand’
     Bitwise AND operator: only bits with values of 1 in both popped
     operands will get the value of 1, the rest will be set to 0.  For
     example, (assuming numbers can be written as bit strings on the
     command-line): ‘00101000 00100010 bitand’ will give ‘00100000’.
     Note that the bitwise operators only work on integer type datasets.

‘bitor’
     Bitwise inclusive OR operator: The bits where at least one of the
     two popped operands has a 1 value get a value of 1, the others 0.
     For example, (assuming numbers can be written as bit strings on the
     command-line): ‘00101000 00100010 bitand’ will give ‘00101010’.
     Note that the bitwise operators only work on integer type datasets.

‘bitxor’
     Bitwise exclusive OR operator: A bit will be 1 if it differs
     between the two popped operands.  For example, (assuming numbers
     can be written as bit strings on the command-line): ‘00101000
     00100010 bitand’ will give ‘00001010’.  Note that the bitwise
     operators only work on integer type datasets.

‘lshift’
     Bitwise left shift operator: shift all the bits of the first
     operand to the left by a number of times given by the second
     operand.  For example, (assuming numbers can be written as bit
     strings on the command-line): ‘00101000 2 lshift’ will give
     ‘10100000’.  This is equivalent to multiplication by 4.  Note that
     the bitwise operators only work on integer type datasets.

‘rshift’
     Bitwise right shift operator: shift all the bits of the first
     operand to the right by a number of times given by the second
     operand.  For example, (assuming numbers can be written as bit
     strings on the command-line): ‘00101000 2 rshift’ will give
     ‘00001010’.  Note that the bitwise operators only work on integer
     type datasets.

‘bitnot’
     Bitwise not (more formally known as one's complement) operator:
     flip all the bits of the popped operand (note that this is the only
     unary, or single operand, bitwise operator).  In other words, any
     bit with a value of ‘0’ is changed to ‘1’ and vice-versa.  For
     example, (assuming numbers can be written as bit strings on the
     command-line): ‘00101000 bitnot’ will give ‘11010111’.  Note that
     the bitwise operators only work on integer type datasets/numbers.


File: gnuastro.info,  Node: Numerical type conversion operators,  Next: Random number generators,  Prev: Bitwise operators,  Up: Arithmetic operators

6.2.4.15 Numerical type conversion operators
............................................

With the operators below you can convert the numerical data type of your
input, see *note Numeric data types::.  Type conversion is particularly
useful when dealing with integers, see *note Integer benefits and
pitfalls::.

   As an example, let's assume that your colleague gives you many single
exposure images for processing, but they have a double-precision
floating point type!  You know that the statistical error a
single-exposure image can never exceed 6 or 7 significant digits, so you
would prefer to archive them as a single-precision floating point and
save space on your computer (a double-precision floating point is also
double the file size!).  You can do this with the ‘float32’ operator
described below.

‘u8’
‘uint8’
     Convert the type of the popped operand to 8-bit unsigned integer
     type (see *note Numeric data types::).  The internal conversion of
     C will be used.

‘i8’
‘int8’
     Convert the type of the popped operand to 8-bit signed integer type
     (see *note Numeric data types::).  The internal conversion of C
     will be used.

‘u16’
‘uint16’
     Convert the type of the popped operand to 16-bit unsigned integer
     type (see *note Numeric data types::).  The internal conversion of
     C will be used.

‘i16’
‘int16’
     Convert the type of the popped operand to 16-bit signed integer
     (see *note Numeric data types::).  The internal conversion of C
     will be used.

‘u32’
‘uint32’
     Convert the type of the popped operand to 32-bit unsigned integer
     type (see *note Numeric data types::).  The internal conversion of
     C will be used.

‘i32’
‘int32’
     Convert the type of the popped operand to 32-bit signed integer
     type (see *note Numeric data types::).  The internal conversion of
     C will be used.

‘u64’
‘uint64’
     Convert the type of the popped operand to 64-bit unsigned integer
     (see *note Numeric data types::).  The internal conversion of C
     will be used.

‘f32’
‘float32’
     Convert the type of the popped operand to 32-bit (single precision)
     floating point (see *note Numeric data types::).  The internal
     conversion of C will be used.  For example, if ‘f64.fits’ is a
     64-bit floating point image, and you want to store it as a 32-bit
     floating point image, you can use the command below (the second
     command is to show that the output file consumes half the storage)

          $ astarithmetic f64.fits float32 --output=f32.fits
          $ ls -lh f64.fits f32.fits

‘f64’
‘float64’
     Convert the type of the popped operand to 64-bit (double precision)
     floating point (see *note Numeric data types::).  The internal
     conversion of C will be used.


File: gnuastro.info,  Node: Random number generators,  Next: Coordinate and border operators,  Prev: Numerical type conversion operators,  Up: Arithmetic operators

6.2.4.16 Random number generators
.................................

When you simulate data (for example, see *note Sufi simulates a
detection::), everything is ideal and there is no noise!  The final step
of the process is to add simulated noise to the data.  The operators in
this section are designed for that purpose.  To learn more about the
definition and implementation "noise", see *note Noise basics::.

   In case each data element's random distribution should have an
independent parameter (for example $\sigma$ in a Gaussian distribution),
the first popped operand can be a dataset of the same size as the
second.  In this case (when the parameter is not a single value, but an
array), each element will have a different parameter.

   When ‘--quiet’ is not given, a statement will be printed on each
invocation of these operators (if there are multiple calls to the
‘mknoise-*’ operators, the statement will be printed multiple times).
It will show the random number generator function and seed that was used
in that invocation.  These are necessary for the future reproducibility
of the outputs using the ‘--envseed’ option, for more, see *note
Generating random numbers::.  For example, with the first command below,
‘image.fits’ will be degraded by a noise of standard deviation 3 units.
     $ astarithmetic image.fits 3 mknoise-sigma

   Alternatively, you can use the operators in this section within the
*note Column arithmetic:: feature of the Table program.  For example,
with the command below, you can generate a random number (centered on 0,
with $\sigma=3$).  With the second command, you can put it into a shell
variable for later usage.

     $ echo 0 | asttable -c'arith $1 3 mknoise-sigma'
     $ value=$(echo 0 | asttable -c'arith $1 3 mknoise-sigma' --quiet)
     $ echo $value

   You can also use the operators here in combination with AWK to easily
generate an arbitrarily large table with random columns.  In the example
below, we will create a two column table with 20 rows.  The first column
will be centered on 5 and $\sigma_1=2$, the second will be centered on
10 and $\sigma_2=3$:

     $ echo 5 10 \
            | awk '{for(i=0;i<20;++i) print $1, $2}' \
            | asttable -c'arith $1 2 mknoise-sigma' \
                       -c'arith $2 3 mknoise-sigma'

   By adding an extra ‘--output=random.fits’, the table will be saved
into a file called ‘random.fits’, and you can change the ‘i<20’ to
‘i<5000’ to have 5000 rows instead.  Of course, if your input table has
different values in the desired column the noisy distribution will be
centered on each input element, but all will have the same
scatter/sigma.

   As mentioned above, you can use the ‘--envseed’ option to pre-define
the random number generator seed (and thus get a reproducible result).
For more on ‘--envseed’, see *note Generating random numbers::.  When
using column arithmetic in Table, it may happen that multiple columns
need random numbers (with any of the ‘mknoise-*’ operators) in one call
of ‘asttable’.  In such cases, the value given to ‘GSL_RNG_SEED’ is
incremented by one on every call to the ‘mknoise-*’ operators.  Without
this increment, when the column values are the same (happens a lot, for
no-noised datasets), the returned values for all columns will be
identical.  But this feature has a side-effect: that if the order of
calling the ‘mknoise-*’ operators changes, the seeds used for each
operator will change(1).

‘mknoise-sigma’
     Add a Gaussian noise with pre-defined $\sigma$ to each element of
     the input dataset (independent of the input pixel value).  $\sigma$
     is the standard deviation of the Gaussian or Normal distribution
     (https://en.wikipedia.org/wiki/Normal_distribution).  This operator
     takes two arguments: the top/first popped operand is the noise
     standard deviation, the next popped operand is the dataset that the
     noise should be added to.

     For example, with the first command below, let's put a Sérsic
     profile with Sérsic index 1 and effective radius 10 pixels,
     truncated at 5 times the effective radius in the center of a mock
     image that is $100\times100$ pixels wide.  We will also give it a
     position angle of 45 degrees and an axis ratio of 0.8, and set it
     to have a total electron count of 10000 (‘1e4’ in the command).
     Note that this example is focused on this operator, for a robust
     simulation, see the tutorial in *note Sufi simulates a detection::.
     With the second command, let's add noise to this image and with the
     third command, we'll subtract the raw image from the noised image.
     Finally, we'll view them both together:

          $ echo "1 50 50 1 10 1 45 0.8 1e4 5" \
                 | astmkprof --mergedsize=100,100 --oversample=1 \
                             --mcolissum --output=raw.fits

          $ astarithmetic raw.fits 2 mknoise-sigma --output=sigma.fits

          $ astarithmetic raw.fits sigma.fits - -g1 \
                          --output=diff-sigma.fits

          $ astscript-fits-view raw.fits sigma.fits diff-sigma.fits

     You see that the final ‘diff-sigma.fits’ distribution was
     independent of the pixel values of the input.  You will also notice
     that within ‘sigma.fits’ the noisy pixels that had a zero value in
     ‘raw.fits’, the noise fluctuates around zero (is negative in half
     of those pixels).  These behaviors will be different in the case
     for ‘mknoise-sigma-from-mean’ below, which is more "realistic" (or
     Poisson-like).

‘mknoise-sigma-from-mean’
     Replace each input element (e.g., pixel in an image) of the input
     with a random value taken from a Gaussian distribution (for pixel
     $i$) with mean $\mu_i$ and standard deviation $\sigma_i$.  Where,
     $\sigma_i=\sqrt{I_i+B_i}$ and $\mu_i=I_i+B_i$ and $I_i$ and $B_i$
     are respectively the values of the input image, and background in
     that same pixel.  In other words, this can be seen as approximating
     a Poisson distribution at high mean values (where the Poisson
     distribution becomes identical to the Gaussian distribution).

     This operator takes two arguments: 1.  the first popped operand
     (just before the operator) is the _per-pixel_ background value (in
     units of electron counts).  2.  The second popped operand is the
     dataset that the noise should be added to.

     To demonstrate the effect of this noise pattern, please run the
     example commands in the description of ‘mknoise-sigma’.  With the
     first command below, let's add this Poisson-like noise (assuming a
     background level of 4 electron counts, to be similar to a
     $\sigma=2$ of the example in ‘mknoise-sigma’).  With the second
     command, let's subtract the raw image from this noise pattern:

          $ astarithmetic raw.fits 4 mknoise-sigma-from-mean \
                          --output=sigma-from-mean.fits

          $ astarithmetic raw.fits sigma-from-mean.fits - -g1 \
                          --output=diff-sigma-from-mean.fits

          $ astscript-fits-view diff-sigma.fits \
                                diff-sigma-from-mean.fits

     You clearly see how the noise in the center of the Sérsic profile
     is much stronger than the outer parts.  As described, above, this
     is behavior we would expect in a "real" observation: the regions
     with stronger signal, also have stronger noise as defined through
     the Poisson distribution
     (https://en.wikipedia.org/wiki/Poisson_distribution)!  The reason
     we described this operator as "Poisson-like" is that, it has some
     shortcomings as opposed to the ‘mknoise-poisson’ operator (that is
     described below):
        • For low mean values (less than 3 for example), this will
          produce a symmetric Gaussian distribution, while the Poisson
          distribution will not be symmetric.
        • The random values from this distribution are floating point
          (unlike the Poisson distribution that produces integers.
        • The random values can be negative (which is not possible in a
          Poisson distribution).

     Therefore to simulate photon-starved images (for example UV or
     X-ray data), the ‘mknoise-poisson’ operator should always be used,
     not this one.  However, in optical (or redder bands) data, the
     background is very bright (much brighter than 10 counts for
     example).  In such cases (as the mean increases), the Poisson
     distributions becomes identical to the Gaussian distribution.
     Furthermore, processed coadded images are no longer integers, but
     floating points with the Sky-level already subtracted (see *note
     Sky value::).  Therefore if you are trying to simulate a processed,
     photon-rich dataset, you can safely use this operator.

     Recall that the background values reported by observatories (for
     example, to define dark or gray nights), or in papers, is usually
     reported in units of magnitudes per arcseconds square.  You need to
     do the conversion to counts per pixel manually.  The conversion of
     magnitudes to counts is described below.  For converting arcseconds
     squared to number of pixels, you can use the ‘--pixelscale’ option
     of *note Fits::.  For example, ‘astfits image.fits --pixelscale’.

     Except for the noise-model, this operator is very similar to
     ‘mknoise-sigma’ and the examples there apply here too.  The main
     difference with ‘mknoise-sigma’ is that in a Poisson distribution
     the scatter/sigma will depend on each element's value.

     For example, let's assume you have made a mock image called
     ‘mock.fits’ with *note MakeProfiles:: and it is assumed zero point
     is 22.5 (for more on the zero point, see *note Brightness flux
     magnitude::).  Let's assume the background level for the Poisson
     noise has a value of 19 magnitudes.  You can first use the
     ‘mag-to-counts’ operator to convert this background magnitude into
     counts, then feed the background value in counts to
     ‘mknoise-sigma-from-mean’ operator:

          $ astarithmetic mock.fits 19 22.5 mag-to-counts \
                          mknoise-sigma-from-mean

     Try changing the background value from 19 to 10 to see the effect!
     Recall that the tutorial *note Sufi simulates a detection:: shows
     how you can use MakeProfiles to build mock images.

‘mknoise-poisson’
     Replace every pixel of the input with a random integer taken from a
     Poisson distribution with the mean value of that input pixel.
     Similar to ‘mknoise-sigma-from-mean’, it takes two operands: 1.
     The first popped operand (just before the operator) is the
     per-pixel background value (in units of electron counts).  2.  The
     second popped operand is the dataset that the noise should be added
     to.

     To demonstrate this noise pattern, let's use ‘mknoise-poisson’ in
     the example of the description of ‘mknoise-sigma-from-mean’ with
     the first command below.  The second command below will show you
     the two images side-by-side, you will notice that the Poisson
     distribution's undetected regions are slightly darker (this is
     because of the skewness of the Poisson distribution).  Finally,
     with the last two commands, you can see the histograms of the two
     distributions:

          $ astarithmetic raw.fits 4 mknoise-poisson \
                          --output=poisson.fits

          $ astscript-fits-view sigma-from-mean.fits poisson.fits

          $ aststatistics sigma-from-mean.fits --lessthan=10
          -------
          Histogram:
           |                       ***
           |                      ******
           |                    **********
           |                    ***********
           |                  **************
           |                 ****************
           |                ******************
           |              **********************
           |             **************************
           |          **********************************  *
           |* **********************************************************
           |------------------------------------------------------------

          $ aststatistics poisson.fits --lessthan=10
          -------
          Histogram:
           |                    *     *
           |                    *     *
           |                    *     *      *
           |             *      *     *      *
           |             *      *     *      *      *
           |             *      *     *      *      *
           |      *      *      *     *      *      *     *
           |      *      *      *     *      *      *     *
           |      *      *      *     *      *      *     *      *
           |      *      *      *     *      *      *     *      *
           |      *      *      *     *      *      *     *      *
           |------------------------------------------------------------

     The extra skewness in the Poisson distribution, and the fact that
     it only returns integers is therefore clear with the commands
     above.  The comparison was further made above in the description of
     ‘mknoise-sigma-from-mean’.  In summary, you should prefer the
     Poisson distribution when you are simulating the following
     scenarios:
        • A photon-starved image (as in UV or X-ray).
        • A raw exposure of a photon-rich image (which may be
          photon-rich, but always integers).

‘mknoise-uniform’
     Add uniform noise to each element of the input dataset.  This
     operator takes two arguments: the top/first popped operand is the
     width of the interval, the second popped operand is the dataset
     that the noise should be added to (each element will be the center
     of the interval).  The returned random values may happen to be the
     minimum interval value, but will never be the maximum.  Except for
     the noise-model, this operator behaves very similar to
     ‘mknoise-sigma’, see the explanation there for more.

     For example, with the command below, a random value will be
     selected between 10 to 14 (centered on 12, which is the only input
     data element, with a total width of 4).

          echo 12 | asttable -c'arith $1 4 mknoise-uniform'

     Similar to the example in ‘mknoise-sigma’, you can pipe the output
     of ‘echo’ to ‘awk’ before passing it to ‘asttable’ to generate a
     full column of uniformly selected values within the same interval.

‘random-from-hist-raw’
     Generate random values from a custom distribution (defined by a
     histogram).  The output will have a double-precision floating point
     type (see *note Numeric data types::).  This operator takes three
     operands:
        • The first popped operand (nearest to the operator) is the
          histogram values.  The histogram is a 1-dimensional dataset (a
          table column) and contains the probability of obtaining a
          certain interval of values.  The histogram does not have to be
          normalized: the GNU Scientific Library (or GSL, which is used
          by Gnuastro for this operator), will normalize it internally.
          The value of each bin (whose probability is given in the
          histogram) is given in the second popped operand.  Therefore
          these two operands have to have the same number of rows.
        • The second popped operand is the bin value (mostly the bin
          center, but it can be anything).  The probability of each bin
          is defined in the histogram operand (first popped operand).
          The bins can have any width (do not have to be evenly spaced),
          and any order.  Just make sure that the same row in the bins
          column corresponds to the same row in the histogram: the
          number of rows in the bins and histogram must be equal.
        • The third popped operand is the dataset that the random values
          should be written over.  Effectively only its size will be
          used by this operator (all values will be over-written as a
          double-precision floating point number).

     The first two operands have to be single-dimensional (a table
     column) and have the same number of rows, but the last popped
     operand can have any number of dimensions.  You can use the
     ‘load-col-’ operator to load the two bins and histogram columns
     from an external file (see *note Loading external columns::).

     For example, in the command below, we first construct a fake
     histogram to represent a $y=x^2$ distribution with AWK. We aim to
     distribute random values from this distribution in a $100\times100$
     image.  Therefore, we use the ‘makenew’ operator to construct an
     empty image of that size, use the ‘load-col-’ operator to load the
     histogram columns into Arithmetic and put the output in
     ‘random.fits’.  Finally we visually inspect ‘random.fits’ with DS9
     and also have a look at its pixel distribution with
     ‘aststatistics’.

          $ echo "" | awk '{for(i=1;i<5;++i) print i, i*i}' \
                          > histogram.txt

          $ cat histogram.txt
          1 1
          2 4
          3 9
          4 16

          $ astarithmetic 100 100 2 makenew \
                          load-col-1-from-histogram.txt \
                          load-col-2-from-histogram.txt \
                          random-from-hist-raw \
                          --output=random.fits

          $ astscript-fits-view random.fits

          $ aststatistics random.fits --asciihist --numasciibins=50
           |                                                 *
           |                                                 *
           |                                                 *
           |                                                 *
           |                                 *               *
           |                                 *               *
           |                                 *               *
           |                *                *               *
           |                *                *               *
           |*               *                *               *
           |*               *                *               *
           |--------------------------------------------------

     As you see, the 10000 pixels in the image only have values 1, 2, 3
     or 4 (which were the values in the bins column of ‘histogram.txt’),
     and the number of times each of these values occurs follows the
     $y=x^2$ distribution.

     Generally, any value given in the bins column will be used for the
     final output values.  For example, in the command below (for
     generating a histogram from an analytical function), we are adding
     the bins by 20 (while keeping the same probability distribution of
     $y=x^2$).  If you re-run the Arithmetic command above after this,
     you will notice that the pixels values are now one of the following
     21, 22, 23 or 24 (instead of 1, 2, 3, or 4).  But the shape of the
     histogram of the resulting random distribution will be unchanged.

          $ echo "" | awk '{for(i=1;i<5;++i) print 20+i, i*i}' \
                          > histogram.txt

     If you do not want the outputs to have exactly the value of the bin
     identifier, but be a randomly selected value from a uniform
     distribution within the bin, you should use ‘random-from-hist’ (see
     below).

     As mentioned above, the output will have a double-precision
     floating point type (see *note Numeric data types::).  Therefore,
     by default each element of the output will consume 8 bytes
     (64-bits) of storage.  This is usually far more than the
     statistical error/precision of your data (and just results in
     wasted storage in your file system, or wasted RAM when a program
     that uses the data is being run, and a slower running time of the
     program).

     It is therefore recommended to use a type-conversion operator after
     this operator to put the output in the smallest type that can be
     used to safely store your data without wasting storage, RAM or
     time.  For the list of type conversion operators, see *note
     Numerical type conversion operators::.  Recall that you already
     know the values returned by this operator (they are one of the
     values in the bins column).

     For example, in the example above, the whole image only has values
     1, 2, 3 or 4.  Since they are always positive and are below 255, we
     can safely place them in an unsigned 8-bit integer (see *note
     Numeric data types::) with the command below (note the ‘uint8’
     after the operator name, and that we are using a different name for
     the output).  After building the new image, let's have a look at
     the sizes of the two images with ‘ls -l’:

          $ astarithmetic 100 100 2 makenew \
                          load-col-1-from-histogram.txt \
                          load-col-2-from-histogram.txt \
                          random-from-hist-raw uint8 \
                          --output=random-u8.fits

          $ ls -lh random.fits random-u8.fits
          -rw-r--r-- 1 name name 85K Jan 01 13:40 random.fits
          -rw-r--r-- 1 name name 17K Jan 01 13:45 random-u8.fits

     As you see, when using a suitable data type, we can shrink the size
     of the file significantly without loosing any information (from 85
     kilobytes to 17 kilobytes).  This difference can be felt much
     better for larger (real-world) datasets, so be sure to always set
     the output data type after calling this operator.

‘random-from-hist’
     Similar to ‘random-from-hist-raw’, but do not return the exact bin
     value, instead return a random value from a uniform distribution
     within each bin.  Therefore the following limitations have to be
     taken into account (compared to ‘random-from-hist-raw’):
        • The number associated with each bin (in the bin column) should
          be its center.
        • The bins have to be in descending order (so the second row in
          the bin column is larger than the first).
        • The bin widths (distance from one bin to another) have to be
          fixed.

     For a demonstration, let's replace ‘random-from-hist-raw’ with
     ‘random-from-hist’ in the example of the description of
     ‘random-from-hist-raw’.  Note how we are manually converting the
     output of this operator into single-precision floating point
     (32-bit, since the default 64-bit precision is statistically
     meaningless in this scenario and we do not want to waste storage,
     memory and running time):
          $ echo "" | awk '{for(i=1;i<5;++i) print i, i*i}' \
                          > histogram.txt

          $ astarithmetic 100 100 2 makenew \
                          load-col-1-from-histogram.txt \
                          load-col-2-from-histogram.txt \
                          random-from-hist float32 \
                          --output=random.fits

          $ aststatistics random.fits --asciihist --numasciibins=50
           |                                          *
           |                                      *** ********
           |                                      ************
           |                                     *************
           |                             *     * *************
           |                         * ***********************
           |                         *************************
           |                         *************************
           |             *************************************
           |********* * **************************************
           |**************************************************
           |--------------------------------------------------

     You can see that the pixels of ‘histogram.fits’ are no longer just
     1, 2, 3 or 4.  Instead, the values within each bin are selected
     from a uniform distribution covering that bin.  This creates the
     step-like feature in the histogram of the output.

     Of course, this extra uniform random number generation can make
     your program slower so be sure to check if it is worth it.  In
     particular, one way to avoid this (and use ‘random-from-hist-raw’
     with a more contiguous-looking output distribution) is to simply
     use a higher-resolution histogram (assuming it is possible: you
     have a sufficient number of data points, or you have an analytical
     expression that you can sample at smaller bin sizes).

     To better demonstrate this operator and its practical usage in
     everyday research, let's look at another example: Assume you want
     to get 100 random star magnitudes that follow the real-world Gaia
     Data release 3 magnitude distribution within a radius of 2 degrees
     around the (RA,Dec) coordinate of (1.23,4.56).  Let's further
     assume that you want to distribute them uniformly over an image of
     size 1000 by 1000 pixels.  So your desired output table should have
     three columns, the first two are pixel positions of each star, and
     the third is the magnitude.

     First, we need to query the Gaia database and ask for all the
     magnitudes in this region of the sky.  We know that Gaia is not
     complete for stars fainter than the 20th magnitude, so we will use
     the ‘--range’ option and only ask for those stars that are brighter
     than magnitude 20.

          $ astquery gaia --dataset=dr3 --center=1.23,3.45 --radius=2 \
                     --column=phot_g_mean_mag --output=gaia.fits \
                     --range=phot_g_mean_mag,-inf,20

     We now have more than 25000 magnitudes in ‘gaia.fits’!  To get a
     more accurate random sampling of our stars, let's construct a
     histogram with 500 bins, and generate our three desired randomly
     selected columns:

          $ aststatistics gaia.fits --histogram --numbins=500 \
                          --output=gaia-hist.fits

          $ asttable gaia-hist.fits -i

          $ echo 1000 \
              | awk '{for(i=0;i<100;++i) print $1/2}' \
              | asttable -c'arith $1 500 mknoise-uniform' \
                         -c'arith $1 500 mknoise-uniform' \
                         -c'arith $1 \
                                  load-col-1-from-gaia-hist.fits-hdu-1 \
                                  load-col-2-from-gaia-hist.fits-hdu-1 \
                                  random-from-hist float32'

     These columns can easily be placed in the format for *note
     MakeProfiles:: to be inserted into an image automatically.

   ---------- Footnotes ----------

   (1) We have defined Task 15971 (https://savannah.gnu.org/task/?15971)
in Gnuastro's project management system to address this.  If you need
this feature please send us an email at ‘bug-gnuastro@gnu.org’ (to
motivate us in its implementation).

