Photometric redshifts for the Pan-STARRS1 survey

We present a robust method to estimate the redshift of galaxies using Pan-STARRS1 photometric data. Our method is an adaptation of the one proposed by Beck et al. (2016) for the SDSS Data Release 12. It uses a training set of 2313724 galaxies for which the spectroscopic redshift is obtained from SDSS, and magnitudes and colours are obtained from the Pan-STARRS1 Data Release 2 survey. The photometric redshift of a galaxy is then estimated by means of a local linear regression in a 5-dimensional magnitude and colour space. Our method achieves an average bias of $\overline{\Delta z_{\rm norm}}=-2.01 \times 10^{-4}$, a standard deviation of $\sigma(\Delta z_{\rm norm})=0.0298$, and an outlier rate of $P_o=4.32\%$ when cross-validating on the training set. Even though the relation between each of the Pan-STARRS1 colours and the spectroscopic redshifts is noisier than for SDSS colours, the results obtained by our method are very close to those yielded by SDSS data. The proposed method has the additional advantage of allowing the estimation of photometric redshifts on a larger portion of the sky ($\sim 3/4$ vs $\sim 1/3$). The training set and the code implementing this method are publicly available at www.testaddress.com.


Introduction
In the last two decades, there has been a rise in the development of large photometric surveys, like the Sloan Digital Sky Survey (SDSS, York et al. 2000), the Panoramic Survey Telescope & Rapid Response System (Pan-STARRS, Chambers et al. 2016), and the Dark Energy Survey (DES, Dark Energy Survey Collaboration et al. 2016). Robust methods to estimate the redshift of galaxies from photometric data are essential to maximize the scientific exploitation of these surveys.
Two main approaches are generally used for the computation of photometric redshifts: methods based on physical models and data-driven methods. In the model-based approach, the estimation of the redshift is obtained by modelling the physical processes that drive the light emission of the object. The simplest and most-commonly used method belonging to this category is spectral energy distribution (SED) fitting. It is based on the definition of a SED model, either from theory or from observations, and the fitting of this model to a series of observations in different bands. The definition of an appropriate model is crucial for the performance of the method, therefore it requires taking many different aspects into account (stellar populations models, nebular emissions, and dust attenuation amongst others). Once the model is defined, observations over the entire wavelength range are required for obtaining an accurate fitting. Examples of these methods are the HYPERZ code (Bolzonella et al. 2000), the BPZ code (Benítez 2000), the LePhare code (Ilbert et al. 2006), and the EAZY code (Brammer et al. 2008). Saglia et al. (2012) also apply a SED technique to compute the photometric redshifts of galaxies using Pan-STARRS broadband photometry.
In large-area photometric surveys like SDSS, Pan-STARRS, and DES, the number of photometric bands available is relatively small (5 for each of the cited surveys) and they only cover the optical part of the spectrum. Thus, if no ancillary data are available, the SED fitting technique is not very robust in the determination of photometric redshifts. On the other hand, these surveys offer a large number of extragalactic sources, being well suitable for the use of data-driven methods. These methods usually employ a supervised machine learning algorithm to estimate the unknown redshift of a galaxy from broadband photometry. Supervised algorithms require a (large) set of reliable spectroscopic redshifts that are used to learn how redshifts correlate with colours. Some examples of these techniques are ANNz (Collister & Lahav 2004), ANNz2 (Sadeh et al. 2016), TPZ (Carrasco Kind & Brunner 2013), GPz (Almosallam et al. 2016), METAPhoR (Cavuoti et al. 2017) or the nearest-neighbors color-matching photometric redshift estimator of Graham et al. (2018).
Another example of a machine learning approach for the computation of photometric redshifts is presented in Beck et al. (2016): a large sample of galaxies (about 2 millions) with both photometric and spectroscopic information is used as training set to estimate the redshift of all the galaxies in SDSS Data Release 12 (DR12, Alam et al. 2015) using a local linear regression. A similar method was also presented in Csabai et al. (2007) and in earlier SDSS releases. Nowadays, these photometric redshifts from SDSS are widely used in a variety of scientific publications and the robustness of the algorithm is well established.
The goal of this paper is to adapt the Beck et al. (2016) algorithm to compute photometric redshifts using the Pan-STARRS1 (PS1) photometric data, which cover an area that is twice as large as the SDSS footprint and whose magnitude limits are about two magnitudes fainter than the SDSS ones. SDSS and PS1 have four photometric bands in common, plus a fifth band that is different, and which is on the bluer side of the spectrum for SDSS and on the redder side for PS1. To adapt Beck et al. (2016) algorithm to PS1 we thus needed to select the appropriate PS1 data that allow us to compute the redshift, to construct a proper training set, and to reassess the performance of the linear regression algorithm when using this information. The training set and the code implementing the PS1 photometric redshift method presented in this paper has been made available for the community at www.testaddress.com.
The method proposed in this paper was initially designed with the purpose of confirming cluster candidates of the Com-PRASS catalogue (Tarrío et al. 2019). This all-sky catalogue of galaxy clusters and cluster candidates was validated by careful cross-identification with previously known clusters, especially in the SDSS and SPT footprints. Still, many candidates remain unconfirmed outside these areas. Having information on the photometric redshifts in the PS1 area will enable us to confirm Com-PRASS candidates in this region. This information will also facilitate the extension of other scientific studies performed with SDSS photometric data to the area of the sky covered by PS1. Some examples are studies related to the formation and evolution of galaxies, or to the properties of dark energy (Salvato et al. 2019).
The paper is organised as follows: Section 2 summarises the linear regression method and how it is adapted to PS1 data. Section 3 describes the procedures that we put in place to prepare the training set using PS1 photometry and SDSS spectroscopy. Section 4 evaluates the performance of the proposed redshift estimation method. Section 5 presents a comparison with the results obtained using different photometric data from PS1 and SDSS. Section 6 gives some practical notes on the use of the method and the associated dataset. Finally, Sect. 7 concludes the paper with a summary of the main results.

Method
In this Section we describe the method to estimate the redshift of galaxies from PS1 photometric data. This method is an adaptation of the method used in the SDSS DR12 (Beck et al. 2016) and can be used to calculate photometric redshifts for all galaxies in the PS1 footprint (∼ 3/4 of the sky). The method is prepared to work on both Data Release 1 (DR1) and Data Release 2 (DR2), although in this paper we will present the results corresponding to DR2. The performance for DR1 data was also tested, finding no significant differences.
The method is data-driven and uses a training set T composed of galaxies with known spectroscopic redshift and a set of magnitudes and colours obtained from the PS1 survey. The redshift of a galaxy is estimated by means of a local linear regression in a D-dimensional magnitude and colour space. The rest of this Section summarizes the linear regression algorithm (see also Beck et al. (2016)) and describes in detail how we selected the magnitude-colour space to be used for PS1. We also explain how we deal with the potential problem of missing information.

Linear regression algorithm
The local linear model establishes that the redshift of a galaxy can be written as a linear combination of D galaxy properties (magnitudes and colours), hereafter features, x 1 , ..., x D , as follows: where x = [1, x 1 , ..., x D ] T is the feature vector of the galaxy. The column vector θ contains the D+1 coefficients of the Ddimensional linear regression, with its first element representing a constant offset. The coefficient vector θ can be esti-mated by constructing an over-determined system of k equations using k galaxies of the training set T : z spec = Xθ, with z spec = [z (1) , ..., z (k) ] T being the spectroscopic redshifts of the k chosen galaxies and X = [x (1) , ..., x (k) ] T the corresponding k feature vectors. The least-squares solution of this system is then: The error of the photometric redshift can be estimated from the difference between the spectroscopic redshifts of the k galaxies and the corresponding photometric redshifts provided by the regression: To apply this method, it is necessary to define how to chose the k training galaxies used to estimate θ, and to define the D features that will characterize each galaxy.
In our case, the k galaxies are chosen to be the nearest neighbours of the target galaxy in terms of Euclidean distance in the D-dimensional space. In particular, we chose k = 100, as in Beck et al. (2016). Additionally, in the case that some of these k neighbours have outlying redshifts (|z ( j) spec −x ( j)Tθ | > 3δ z phot ), we discard them and repeat the computation ofθ (Eq. 2) using the remaining l < k neighbours. We note that, in some cases, a galaxy can fall outside the D-dimensional bounding box of its nearest neighbours. In these cases, Eq. 1 constitutes an extrapolation, so the results may be less reliable. The impact of the extrapolation on the estimated photometric redshift is evaluated in Section 4.2. Our code provides a flag to indicate these cases.

Feature selection
The key point to successfully employ the linear regression algorithm described above with PS1 photometric data is to appropriately select the D features to be used. SDSS and PS1 surveys have both imaged the sky using five broadband filters. Four of these filters (g, r, i, and z) are similar in both surveys, although with some minor differences . The fifth filter, however, is completely different: SDSS uses the u filter, which covers the bluest part of the measured spectrum (at bluer waveleghts than the g filter), whereas PS1 uses the y filter, which spans the reddest part of the spectrum (at redder wavelengths than the z filter). The method defined in Beck et al. (2016) uses the SDSS r magnitude, and the u − g, g − r, r − i, and i − z colours to define the 5-dimensional space in which the linear regression takes place to estimate the redshift. Since the u-band is not available in PS1, a natural choice inspired in Beck et al. (2016) is to use the PS1 r magnitude, and the four colours that can be constructed with consecutive magnitudes, i.e., g − r, r − i, i − z, and z − y. These are the 5 features that we decided to use in our method. However, it is worth noticing that other combinations of the 5 bands are also possible without significant difference in the results, given that all the photometric information is included.
PS1 database provides several ways of measuring magnitudes and fluxes of objects in its five photometric bands. We used stack photometry since it provides the best signal-to-noise, according to Magnier et al. (2019). Then, different photometric measurements are available: -PSF magnitudes: obtained from fitting a predefined PSF form to the detection. These magnitudes are especially relevant for point sources (e.g. stars).
-Kron magnitudes: inferred from the growth curve, after determining the Kron radius of the object. These magnitudes are especially relevant for non-point-sources. -Aperture magnitudes: they measure the total count rate for a point source based on integration over an aperture plus an extrapolation involving the PSF. According to the PS1 database documentation, this photometry should not be used for extended sources, so we will not use it in our method. -Fixed-aperture measurements: flux measured within several predefined aperture radii (1.03, 1.76, 3.00, 4.63, and 7.43 arcsec).
Kron magnitudes are the most appropriate for extended objects like galaxies, so we have chosen to use the r-band Kron magnitude for defining our r feature. Regarding the four colour features (g − r, r − i, i − z, z − y), we have considered two different approaches: they can be calculated either a) from fixed-aperture fluxes, or b) from Kron magnitudes.
The first approach (aperture colours) computes the four colours within a fixed aperture. To obtain the aperture magnitudes within the most appropriate aperture, we selected for each galaxy the g, r, i, z, and y fixed-aperture fluxes corresponding to the closest aperture to the r-band Kron radius of the galaxy (rkronrad). Then, the five selected aperture fluxes are converted into aperture magnitudes.
The full PS1 dataset files available for direct download do not provide the above-mentioned fixed-aperture fluxes, which need to be queried to the database. Instead, they provide Kron magnitudes. As an alternative approach, we evaluated the use of these magnitudes to compute the four colours required by our method. We note that the colours constructed from the Kron magnitudes are not physically motivated, since the five different magnitudes are not measured within the same aperture. However, we will show that they provide very similar results to the ones obtained when using the fixed-aperture colours defined above, so for convenience, we added this alternative in our code. Unless otherwise stated, the results presented in this paper were obtained with the aperture colours calculated from the fixed-aperture fluxes. We include a comparison between the different approaches in Sect. 5.

Feature computation
Before calculating the five features we need to apply a dereddening correction to the downloaded or calculated magnitudes. Reddening is produced by the scattering of the light by dust in the interstellar medium and it depends on the position of the object in the sky. Therefore, it has to be corrected in order to obtain magnitudes that are more correlated with redshift.
We obtained this correction in the following way: Firstly, we computed the colour excess E(B-V) for each galaxy using the Schlegel et al. (1998) maps. Then, we obtained the extinction A λ for the g, r, i, and z bands by multiplying the colour excess by the values presented in table 22 of Stoughton et al. (2002) for the g, r, i, and z SDSS filters respectively, which are very similar to the ones used in PS1. For the y band, which is not present in SDSS, we calculated the extinction using the parametrization of Fitzpatrick (1999) taking the effective λ of the y band (λ eff = 9620 Å) presented in table 4 of Tonry et al. (2012).
We then applied the dereddening correction (g = g downloaded − A g , and equivalently for the other bands), and we computed the five features (g − r, r − i, i − z, z − y, and r).
Each dimension is then standardized, by removing the mean and dividing by the standard deviation of the training set T . In this way, all the features span similar ranges and, thus, contribute with a similar weight to the linear regression. This feature scaling is a common practice in algorithms which use Euclidean distance, like ours, since otherwise the feature with a larger scale (the magnitude in our case) would dominate the computation of the distance.
We note that the zero-point correction that is usually taken into account in public software for the computation of photometric redshifts does not need to be included in our method. The reason is that our method is not sensitive to the addition of constant terms to the features, since it uses standardized features.

Missing features
The PS1 dataset contains galaxies for which one or more magnitudes may not be available, resulting in missing features. The method described above can still be applied to estimate the photometric redshift of these galaxies in several ways. In particular, we have decided to calculate the redshift of such galaxies by using only the available features, i.e. we construct the feature vector x with D < D features, both for the target galaxy and the training galaxies, and then use Eqs. 1 and 2 as before. In this way, we will use the subset of the training set T that has all the five features available (T 5 ) to calculate the redshift of a galaxy that has the five features. Likewise, when a galaxy is missing one feature, we will also use T 5 as training set, but we will not consider the missing feature in any of the training galaxies. Another possible approach to follow in the case of a missing feature in the target galaxy could be to use as training set the subset of T that has the other 4 features available (T 4 , with T 5 ⊂ T 4 ⊂ T ). In this paper, we report the results corresponding to the first approach, but we note that the second option produces similar results and is also available in the code.
It is worth mentioning that, for simplicity, the standardization of the features is done in any case with the mean and the standard deviation of the subset T 5 . We tested that other reasonable choices (e.g. using, for each feature, the mean and the standard deviation of all the galaxies containing that feature) do not yield any significant difference in the results.

Training set
The training set of Beck et al. (2016) included more than 2 millions galaxies with spectroscopic redshifts. In this work, our goal was to use the same training set, but with features obtained from PS1 magnitudes instead of SDSS. To construct it, we made use of the CasJobs tool in SDSS, which allows querying both the SDSS and the PS1 databases. In the catalogue of Beck et al. (2016) each galaxy is identified via its ObjID (a unique number assigned to each object in the SDSS database). Thus, our first step was to obtain the coordinates for each object using an appropriate query to the SDSS database. Then, we performed a query in the PS1 database  to look for the PS1 object nearest to each SDSS object. This was done using the fGetNearestObjEq function, and limiting the search to a radius of 30 . This conservative choice allowed us to define a posteriori the matching distance up to which we can consider the match reliable. Objects with greater matching distances are not kept in the training set. We will show later that this maximum distance was fixed to 1 .
Our query produces the following output parameters: the identifier (ObjID) and coordinates (ra, dec) of the object in both the SDSS and PS1 databases, the distance between the two positions, the g, r, i, z, and y Kron magnitudes and their associated errors in the PS1 database (using stack photometry), the PS1 primarydetection flag, the Kron radius measured in the PS1 r band, and the g, r, i, z, and y PS1 fluxes measured within the five predefined aperture radii (1.03, 1.76, 3.00, 4.63, and 7.43 arcsec) and their associated errors. Table 1 lists these parameters and the tables where they are available. The training set T is constructed from this catalogue, after cleaning it for unwanted objects, calculating the five features defined in Sect. 2, and taking the spectroscopic redshift from the catalogue of Beck et al. (2016). In the following Subsections, we describe these steps in detail.

Cleaning
The catalogue resulting from the query to the PS1 database contains some duplicate entries, i.e. objects with the same ObjID and the same or different properties. We cleaned this catalogue from these objects by keeping only one object for each ObjID. In particular, we selected the ones for which primarydetection was equal to 1, which indicates that the entry is the primary stack detection. If there was more than one object satisfying this condi- tion, we selected the one with more magnitudes available. Exact duplicates were also removed.
We additionally removed from the downloaded catalogue two classes of objects. The first class corresponds to objects for which PS1 and SDSS photometry are very different. PS1 and SDSS have four photometric bands in common (g, r, i, and z), so we expect to have a small difference between the magnitudes in those four bands measured by PS1 and SDSS. However, we noticed that our catalogue included some objects for which the difference between these magnitudes was very high (even more than ten magnitudes in some cases). For precaution, we decided to exclude them from our training set. In particular, we excluded objects for which the difference between any SDSS magnitude and the corresponding PS1 magnitude is greater than 4.
The second class of excluded objects corresponds to those that appear to be too bright for their assigned spectroscopic redshift. We noticed the presence of very bright objects at high redshift in our catalogue that are not physically possible (e.g. r = 12.1 at z = 0.82). After a visual inspection in SDSS we found that these objects were indeed bad samples due to two main reasons: a) low-redshift star-forming galaxies with wrong (high) SDSS spectroscopic redshift; and b) wrong magnitude measurements (for example, in galaxies affected by the light of close-by saturated stars, or by external regions of foreground ex- Table 2. Magnitude limits for different redshift ranges. Objects with magnitudes below these limits are not included in the training set.

Magnitude limit
Redshift range g = 11.00 + 20.00 z spec 0.1 < z spec < 0.2 g = 12.33 + 13.33 z spec 0.2 < z spec < 0.5 g = 17.75 + 2.50 z spec 0.5 < z spec < 0.9 g = 20.00 z spec > 0.9 r = 10.00 + 20.00 z spec 0.1 < z spec < 0.2 r = 12.00 + 10.00 z spec 0.2 < z spec < 0.5 r = 14.00 + 6.00 z spec 0.5 < z spec < 1.0 r = 20.00 z spec > 1.0 i = 10.00 + 20.00 z spec 0.1 < z spec < 0.2 i = 12.00 + 10.00 z spec 0.2 < z spec < 0.4 i = 14.00 + 5.00 z spec 0.4 < z spec < 1.0 i = 19.00 z spec > 1.0 z = 10.75 + 12.50 z spec 0.1 < z spec < 0.3 z = 12.25 + 7.50 z spec 0.3 < z spec < 0.5 z = 14.00 + 4.00 z spec 0.5 < z spec < 1.0 z = 18.00 z spec > 1.0 y = 9.25 + 17.50 z spec 0.1 < z spec < 0.3 y = 12.25 + 7.50 z spec 0.3 < z spec < 0.5 y = 14.00 + 4.00 z spec 0.5 < z spec < 1.0 y = 18.00 z spec > 1.0 tended galaxies). We decided to remove these objects by setting, for each magnitude, the limits in the magnitude-redshift relations given in Table 2. Fig. 1 shows the magnitude limits for the rband, together with the galaxies that are removed after applying the different magnitude limits. Figure 2 shows the distribution of distances between the PS1 and SDSS objects in the cleaned catalogue. The scale is logarithmic to highlight that there is a small tail of objects far away from the SDSS position, as allowed from our large search radius (30 ). Since our goal is to keep only the objects for which the match is secure, we removed from our sample all the objects whose distances to the SDSS positions were larger than 1 . These were 2368 objects out of 2 316 092 (corresponding to ∼ 0.10%), resulting in a training set T with 2 313 724 objects. This conservative approach allows us to safely use the SDSS spectroscopic redshift with the PS1 magnitudes.

Final training set
The final training set T contains 2 313 724 galaxies. For each galaxy, the training set provides the spectroscopic redshift z spec obtained from the catalogue of Beck et al. (2016), and the 5 features (g − r, r − i, i − z, z − y, and r) obtained as explained in Sect. 2.3. The features that could not be calculated due to a missing magnitude were set to a default value (-999).
The redshift distribution of the galaxies in T is shown in Fig. 3, where T has been divided into two subsets: the galaxies that have the 5 features available (T 5 , in blue), and the galaxies that have one or more features missing (in red). For comparison, Fig. 3 also shows the spectroscopic redshift distribution of the original SDSS training set, which contained 2 379 096 galaxies.
The PS1 training set that we are presenting in this Section is smaller than the original SDSS one, as expected when doing a match between different catalogues. This difference can be due to errors in the astrometry or photometry of the two surveys, as Fig. 3. Distribution of the spectroscopic redshifts in our training set. In blue we show the galaxies that have the 5 features available (T 5 ). In red, the remaining galaxies. The dashed black line shows the redshift distribution of the original SDSS training set. well as to intrinsic limits of our match methodology. However, we stress that we tried to use a conservative approach in which the number of galaxies is smaller, but the robustness of the match is favoured. This result was reached while loosing less than 3% of the original galaxies.

Overall redshift precision
To evaluate the performance of the proposed method, we randomly divided the training set T into two disjoint subsets of different sizes: T train and T test . The former is used as training set for estimating the photometric redshift of the galaxies in the latter. We chose the sizes to be 99% and 1% of T , respectively. We repeated the experiment 10 times to analyse the stability of the results. Figure 4 shows the photometric redshift z phot , the actual error z phot − z spec , and the error divided by the estimation of the error provided by the method (z phot − z spec )/δ z phot as a function of the spectroscopic redshift for the galaxies in the 10 different T test sets having the 5 magnitudes available. The photometric redshift follows quite well the spectroscopic redshift, especially in the intermediate redshift range (0.1 < z spec < 0.6), where the average bias ∆z = |z phot − z spec | is below 0.02 or 0.5δ z phot .
For high-redshift galaxies (z spec > 0.6), the method tends to underestimate the redshift, and also presents a higher scatter. This behaviour was also observed in Beck et al. (2016). The increased scatter is due to the low number of high-redshift galaxies in the training set. The negative bias is an Eddington bias produced by the limited depth of the PS1 survey: close to the detection limit, over-luminous galaxies are preferentially detected, yielding a bias towards lower redshifts. In this redshift range 54% of the galaxies are within ±δ z phot and 86% are within ±2δ z phot . The percentage of galaxies in this range whose redshift is estimated via extrapolation is 0.10%, six times higher than the value in the intermediate redshift range (0.017%). For lowredshift galaxies (z spec < 0.1), the method tends to overestimate the redshift, as in Beck et al. (2016). In this redshift range 65% of the galaxies are within ±δ z phot and 94% are within ±2δ z phot . The percentage of galaxies in this range whose redshift is estimated via extrapolation is 0.074%, higher than in the intermediate redshift range, but lower than in the high redshift range. Fig. 4. Photometric redshift z phot (left), redshift error z phot − z spec (middle), and error divided by the estimation of the error provided by the method (z phot − z spec )/δ z phot (right) as a function of the spectroscopic redshift z spec . The black dots represent each individual galaxy of the 10 different T test sets. Only the galaxies with the 5 magnitudes available were included. The red solid and dotted lines represent the median and the 68% confidence regions, respectively, computed for groups of 1000 galaxies with consecutive z spec . The orange line shows z phot = z spec . In order to quantitatively compare the average performance of the proposed method to the one obtained in Beck et al. (2016), we use the same definition of the normalized redshift estimation error, that is ∆z norm = z phot −z spec 1+z spec . After iteratively removing the outliers, defined as |∆z norm | > 3σ(∆z norm ), the average bias of our method is ∆z norm = −2.01 × 10 −4 , the standard deviation is σ(∆z norm ) = 0.0298, and the outlier rate is P o = 4.32%, when calculated on the ensemble of results from the 10 experiments. The differences between the 10 experiments were negligible. For reference, the results reported by Beck et al. (2016) were ∆z norm = 5.84 × 10 −5 , σ(∆z norm ) = 0.0205, and P o = 4.11%, which are of the same order, but slightly better than ours. We note however that they were calculated using a different training set, so they are not directly comparable. A direct comparison is presented in Sect. 5.2. Figure 5 shows the normalized histogram of (z phot − z spec )/δ z phot together with a standard normal distribution. The two distributions are well in agreement, apart from a small bias (as in Beck et al. 2016, cfr. their fig. 4). This indicates that the estimated errors δ z phot represent quite well the accuracy of the redshift estimation.

Impact of the photometric errors
The errors in the measurements of the photometric magnitudes have an impact in the final accuracy of the estimated redshift. To evaluate this effect, we classified the galaxies into five different classes according to their photometric errors. Class 1 includes galaxies with low photometric errors, and classes 2-5 include galaxies with progressively higher errors. The error limits for the different classes were manually chosen and are given in Table 3. We also define an additional class E, which includes the galaxies whose redshift is estimated via an extrapolation of Eq. 1. This occurs when the galaxy features lie outside the bounding box of its nearest neighbours, as mentioned in Sect. 2.1.
The photometric error corresponding to the r-band Kron magnitude (∆r) is directly obtained from the query to the PS1 database (rKronMagErr in Table 1). The photometric errors in the four aperture colours are obtained from the errors in the corresponding aperture fluxes. If f g is the aperture flux in the g band and ∆ f g the corresponding error, which are obtained from the query (rgc6flxR and rgc6flxErrR in Table 1), the error in the g-band aperture magnitude is calculated as ∆g = 2.5 log(e) × ∆ f g / f g , and analogously for the other aperture magnitudes. The error in the aperture colours is thus calculated as ∆(g − r) = (∆g) 2 + (∆r) 2 , and similarly for the other colours. Table 4 summarizes the performance of the redshift estimation for the different photometric classes. The bias is very close to 0 for all the classes, being positive for class 1 and increasingly negative for classes 2 to 5. This results in a slightly negative bias for the whole sample (∆z norm = −2.01 × 10 −4 ) since, even though the number of class 1 galaxies dominates, the negative bias of the other classes is higher in absolute value. The standard deviation of the normalized redshift estimation error σ(∆z norm ) increases for higher photometric errors, as expected, and the same occurs with the outlier rate. For the galaxies in class E, the bias and σ(∆z norm ) is higher than for the other classes.
The defined photometric classes can be used to filter out, if needed, the galaxies for which the redshift estimation is less precise.

Impact of the position in the colour-magnitude space
The position of the galaxy in the D-dimensional feature space also has an effect on the redshift estimation error. Galaxies situated in dense regions are expected to have smaller errors, since their neighbours will be very close to them in the D-dimensional colour-magnitude space, and will probably have similar red-P. Tarrío and S. Zarattini : Photometric redshifts for the Pan-STARRS1 survey Fig. 6. Photometric redshift results in the 10 different T test sets as a function of the r − i and the g − r colours. The left panel shows the average standard deviation of the redshifts of the nearest neighbours σ(z NN ), the middle panel shows the rms of the actual error z phot − z spec , and the right panel shows the average estimated errors δ z phot . For easier comparison, the scale in the three panels was set between 0 and 0.1, with the red colour indicating errors that are bigger or equal than 0.1. For reference, the black lines represent the contours of the galaxy count distribution of the training set T 5 , with the four displayed contours corresponding to 1000, 300, 100 and 10 galaxies per colour bin.  Table 3. Photometric error limits for the defined classes. A galaxy belongs to a given class if its five photometric errors are below the specified values for that class and it is the lowest possible class. Class 5 contains galaxies for which one or more of the photometric errors are above the limits corresponding to class 4. shifts. On the contrary, galaxies situated in sparse regions will have larger errors in the redshift estimation, since their neighbours will be further away in the colour-magnitude space, and will probably have a bigger dispersion in their redshifts.
To characterize this effect we have computed several error maps that provide the redshift estimation errors as a function of the position in the D-dimensional colour-magnitude space. These error maps can be used to filter out, if required, the regions in the colour-magnitude space that have larger errors. Figure 6 illustrates this effect in the g−r and r−i colour plane. The colour maps in this figure show three measurements of the redshift estimation error as a function of g − r and r − i: the average standard deviation of the redshifts of the nearest neighbours σ(z NN ), the root mean square (rms) of the actual error z phot −z spec , Table 4. Average normalized redshift estimation bias ∆z norm , standard deviation σ(∆z norm ) and outlier rate P o for the defined photometric classes. These quantities were calculated after iteratively removing the outliers, defined as |∆z norm | > 3σ(∆z norm ). The number of galaxies N, from the 10 different T test sets, belonging to each class is also indicated. and the average estimated errors δ z phot . The different error measurements show a similar behaviour. The estimated error δ z phot is closely related to the actual error, which further supports that it is a good estimator of the error, as previously shown in Fig. 5. As expected, both errors are clearly correlated with the deviation of the redshifts of the nearest neighbours. In the regions where the dispersion is higher, the redshift estimation has a bigger error. The contour lines in Fig. 6 represent the galaxy count distribution of the training set T 5 . By comparing these contours with the Article number, page 7 of 12 A&A proofs: manuscript no. Paper-tarrio  Fig. 8. Left panel: Average normalized error ∆z norm = (z phot − z spec )/(1 + z spec ) in the 10 different T test sets as a function of the magnitude r and the spectroscopic redshift z spec . Right panel: redshift error z phot − z spec as a function of the spectroscopic redshift z spec for galaxies with 18 < r < 20 (orange dots) and 20 < r < 21 (black dots). Each dot represents an individual galaxy of the 10 different T test sets. The thick solid and dotted lines represent the median and the 68% confidence regions, respectively, computed in small z spec intervals for the galaxies with 18 < r < 20 (blue lines) and 20 < r < 21 (red lines). The green line shows z phot = z spec .
background error maps, we see that there is a clear correlation between the photometric redshift errors and the galaxy count distribution: denser regions yield smaller errors and sparser regions yield bigger errors. Figure 7 shows the same three measurements of the redshift estimation error shown in Fig. 6 as a function of r and i − z. The behaviour is similar to that observed in Fig. 6, with an estimated error δ z phot closely following the deviation of the redshifts of the nearest neighbours σ(z NN ) and the actual error z phot − z spec . The contour lines in this figure show the galaxy count distribution as a function of r and i − z. By comparing these contours with the error maps we see again a correlation between the two, although less clear than in Fig. 6. This is due to an additional effect that is analysed next: fainter galaxies tend to have larger errors than brighter galaxies.
In fact, one of the input features of the proposed method depends directly on the r-band magnitude of the galaxies, and thus, on their apparent brightness. This feature has a strong impact on the estimation of the photometric redshift, as shown in Fig.  8. The left panel of Fig. 8 shows the average normalized error ∆z norm = (z phot − z spec )/(1 + z spec ) as a function of the magnitude r and the spectroscopic redshift z spec . While for brighter galaxies (r < 20) the average normalized error is below 0.1, for fainter galaxies (r > 20) the error increases significantly, especially for galaxies at z < 0.4 or z > 0.8. The right panel of Fig. 8 shows the redshift error z phot − z spec for bright (18 < r < 20) and faint (20 < r < 21) galaxies. While the redshift of brighter galaxies is well estimated, with a small bias both at low and high redshift, the error for fainter galaxies is higher, especially when the true redshift is far from 0.5 − 0.6.

Impact of missing features
The proposed method is also able to work when one or several features are missing. When this occurs, the performance of the method degrades, with an increased scatter in the photometric versus spectroscopic redshift relation. Missing features usually appear in very faint galaxies, but can also occur in brighter galaxies due to photometric measurement errors. In our training set T with 2 313 724 galaxies, the r Kron magnitude is missing from Table 5. Average normalized redshift estimation bias ∆z norm , standard deviation σ(∆z norm ) and outlier rate P o for the experiments in which one of the features is removed. These quantities were calculated after iteratively removing the outliers, defined as |∆z norm | > 3σ(∆z norm ). 48 416 galaxies (2.1%), and the g − r, r − i, i − z and z − y aperture colours are missing from 144 352 (6.2%), 51 150 (2.2%), 52 357 (2.3%), and 57 350 (2.5%) galaxies, respectively. Most of these galaxies are faint. For example, approximately 91% of the galaxies without the g − r aperture colour have an r Kron magnitude r > 20, while only 9% are brighter than r = 20.
To evaluate the effect of a missing feature independently of the position of the galaxy in the magnitude-colour space, we artificially removed one of the features from our training set T 5 , and repeated the experiment described in Section 4.1, using the four remaining features for both the test and training subsets. The results are similar to those presented in Fig. 4, but with a higher scatter. Table 5 summarizes the results obtained when the different features are removed. When the r Kron magnitude is removed, the standard deviation of the normalized bias is σ(∆z norm ) = 0.0365, 22% higher than when using the 5 features (σ(∆z norm ) = 0.0298). The effect is smaller when one of the aperture colours is removed, with an increase of 13%, 11%, 4%, and 1% in σ(∆z norm ) for g − r, r − i, i − z, and z − y, respectively. This indicates that, among the five features, the r Kron magnitude has the strongest effect in the determination of the photometric redshifts, while the aperture colours play a weaker role. On the other hand, the average bias remains small, and the outlier rate is not much affected by the removal of any of the colour features, but increases when the r magnitude is removed.  Fig. 9. Average normalized error ∆z norm = (z phot −z spec )/(1+z spec ) in the 10 different T test sets as a function of the magnitude r and the spectroscopic redshift z spec , for three different sets of features: left panel corresponds to PS1 features with aperture colours; middle panel corresponds to PS1 features with Kron colours; and right panel corresponds to SDSS features. In the three cases the training set is T Beck 9 .

Comparison with other photometric features
The photometric redshift estimation method described in Sect. 2 is a general technique that can be applied to different sets of features. In the previous Section we presented the results obtained when using the five PS1 features described in Sect. 2.2, i.e., the PS1 r-band Kron magnitude and the g − r, r − i, i − z, and z − y aperture colours. In this Section we analyse the effects of using different sets of features. In particular, we consider two different cases: 1) PS1 Kron colours, and 2) SDSS features, as in Beck et al. (2016). In the first case, we will assess whether the Kron colors, which are not physically motivated (see Sect. 2.2) but directly available for download for the complete PS1 survey, can be used if more convenient. Considering the SDSS features will allow us to compare the performance of the method using PS1 information with respect to the original method of Beck et al. (2016). In order to do a fair comparison, we restricted our training set to the galaxies in T that have the r-band Kron magnitude, the four aperture colours, and the four Kron colours available in the PS1 dataset. Moreover, we also discard the galaxies that do not satisfy the colour cut and photometric error criteria defined in Eq. 7 of Beck et al. (2016), based on SDSS information. The resulting training set, T Beck 9 , contains 1 776 508 galaxies. As in the previous experiments described in Sect. 4, we divided 10 times T Beck 9 into a training subset and a test subset. Then, we computed the photometric redshift of the galaxies in the 10 test sets using as features: 1) the r-band Kron magnitude and the four aperture colours from PS1; 2) the r-band Kron magnitude and the four Kron colours from PS1; and 3) the five SDSS features defined in Beck et al. (2016). Table 6 reports the average bias, standard deviation, and outlier rate obtained in the three cases. Figure 9 shows the average normalized bias in the three cases as a function of the r-band Kron magnitude and the spectroscopic redshift z spec . The results obtained with PS1 aperture colours are not exactly the same as in the experiment presented in Sect. 4 (see Table 4 and Fig. 8) because the training set (T Beck 9 instead of T 5 ) contains now less galaxies. The galaxies that have been removed with respect to T 5 are those for which a PS1 Kron magnitude is missing or that do not satisfy the SDSS criteria defined in Beck et al. (2016). These are probably galaxies with poorer photometry, which explains the slight improvement in the performance (lower standard deviation and outlier rate) with respect to the results presented in Sect. 4.

Aperture colours vs Kron colours
As shown in Table 6 and Fig. 9, the results obtained when using Kron colours are very similar to the ones obtained when using aperture colours. We have also seen a nearly identical performance to the one shown in Figs. 4 and 5: the photometric redshift calculated from Kron colours follows quite well the spectroscopic redshift, especially in the intermediate redshift range (0.1 < z spec < 0.6), and the method tends to underestimate the redshift for high redshift galaxies (z spec > 0.6) and to overestimate the redshift for low redshift galaxies (z spec < 0.1).
Since the difference between using PS1 aperture colours or PS1 Kron colours is negligible, we have included the possibility of selecting which features to use in the code available at www. testaddress.com. The user may choose the one that is more convenient, without significant impact on the results. Figure 9 and Table 6 show that using SDSS information instead of PS1 information results in a slightly better performance. The standard deviation of the normalized redshift error is lower with SDSS features. Moreover, although the global average bias is smaller for PS1 (with aperture colours), Fig. 9 shows that SDSS features result in a lower bias both at high and low redshift. This effect is especially noticeable for fainter galaxies. In the following, we analyse the possible causes of this behaviour.

PS1 features vs SDSS features
Firstly, PS1 and SDSS features are defined differently. One of the SDSS features is the u − g colour, which is not available in Pan-STARSS; whereas the z−y colour is available in PS1 but not in SDSS. Including a bluer information in SDSS allows a better estimation of the redshift of lower redshift galaxies. To check if this is enough to explain the observed behaviour, we repeated the calculation of the photometric redshift removing the u − g colour Article number, page 9 of 12 A&A proofs: manuscript no. Paper-tarrio from SDSS features and removing the z − y colour from PS1 features. The results with SDSS features degrade, but still show a slightly better performance than with PS1 features, so this feature difference does not entirely explain the better behaviour of SDSS features.
Secondly, the magnitudes g, r, i and z have different values in SDSS and PS1. It turns out that SDSS magnitudes show a better correlation with the spectroscopic redshift, which explains their power to better estimate the photometric redshift. Figure 10 shows a comparison of the correlation between the r − i colour and the spectroscopic redshift for SDSS and PS1 (considering the r − i aperture colour). For a given redshift, PS1 values show a larger dispersion than SDSS. The same occurs for the g − r and i − z aperture colours, for the Kron colours, and for the r Kron magnitude. The reason for this larger scatter could be that PS1 aperture colours are not measured exactly at the Kron radius, but at the closest one from the five available apertures, resulting in colours that do not correspond to the same percentage of flux for all the galaxies. Conversely, SDSS colours are computed from ModelMag magnitudes, so they correspond to the total flux of the galaxy. On the other hand, PS1 Kron colours are not physically motivated, since they are calculated as the difference between two magnitudes that may be measured in different radii.

Practical guidelines for using the method
The training set T and the code implementing our method are available for download at the following webpage: www. testaddress.com. This allows the estimation of the photometric redshift of any galaxy in the PS1 survey. The code includes several configuration options that are described in detail in the webpage. The two main options are the choice between aperture or Kron magnitudes, and the choice of the subset of T to be used for training.
Depending on the required accuracy and on the specific use of the photometric redshifts provided by our method, one may want to use all the possible photometric redshifts regardless of their error, or prefer to use a lower amount of more accurate photometric redshifts. In this Section, we summarize the different ways that we provide of selecting the best redshifts.
There are three main parameters that can be used for this selection: the estimated error δ z phot , the photometric error class, and the extrapolation flag. Figure 11 shows the effect of using different cuts in the photometric redshift errors. As expected, introducing a cut in δ z phot reduces the errors (see for comparison Fig. 4, where no cuts were used). However, if the cut is too severe, the resulting sample may be limited in terms of redshift and colour space coverage. Therefore, we suggest to test different values to find the most appropriate for a particular goal. Figure 12 shows the effect of adding a cut using the photometric error class. By comparing Fig. 11 and Fig. 12 we can see that the photometric error class selection mainly reduces the scatter at high redshifts, where the photometric errors are larger.
Given the small number of galaxies in our dataset for which an extrapolation was performed, filtering them out does not bring a noticeable effect on the ensemble. However, this parameter is a good indicator of the accuracy of the results, as shown in Table  4, so it can be used to filter out some of the calculated redshifts when there is a need for high accuracy.
Finally, the error maps presented in Sect. 4.3 can be used to filter out some of the galaxies located in the regions of the magnitude-colour space that are more prone to errors.

Summary
We present a data-driven method to compute photometric redshifts for galaxies using the PS1 survey. In this work we use data from the PS1 DR2, but we tested the results also for the DR1, finding no significant difference. Our method is an adaptation of the one proposed by Beck et al. (2016) for the SDSS DR12, based on a local linear regression in a 5-dimensional magnitude and colour space. To adapt Beck et al. (2016) algorithm to PS1 we select appropriate magnitudes and colours (r, g − r, r − i, i − z, and z − y) for defining the 5-dimensional space, and we construct a proper and clean training set composed of 2 313 724 galaxies, whose spectroscopic redshift is available from SDSS and whose magnitudes and colours are obtained from the PS1 DR2 survey.
We assess the performance of this method by means of a cross-validation on the training set, i.e., we use part of the galaxies of our training set as test galaxies to estimate their photometric redshifts and we then compare them to their true (spectroscopic) redshifts. We estimate that the average bias of Fig. 11. Photometric redshift as a function of spectroscopic redshift, for three different sets. The red solid and dashed lines represent the median and the 68% confidence regions, respectively, computed in small z spec intervals. The orange line shows z phot = z spec . On panel (a), we included the galaxies with a reported redshift error of δ z phot < 0.05. On panel (b), we included the galaxies with a reported redshift error of δ z phot < 0.03. On panel (c), we included the galaxies with a reported redshift error of δ z phot < 0.02. Fig. 12. Photometric redshift as a function of spectroscopic redshift, for three different sets. The red solid and dashed lines represent the median and the 68% confidence regions, respectively, computed in small z spec intervals. The orange line shows z phot = z spec . On panel (a), we included the galaxies in photometric error class 1 and with a reported redshift error of δ z phot < 0.05. On panel (b), we included the galaxies in photometric error class 1 and with a reported redshift error of δ z phot < 0.03. On panel (c), we included the galaxies in photometric error class 1 and with a reported redshift error of δ z phot < 0.02. our method is ∆z norm = −2.01 × 10 −4 , its standard deviation, σ(∆z norm ) = 0.0298, and the outlier rate P o = 4.32%.
We also evaluate the impact of the photometric uncertainties on our redshift determination. This was done by dividing the entire sample in 5 photometric classes of growing photometric errors. As expected, the uncertainties on the photometric redshifts are smaller where the photometric errors are smaller. There is also a fraction of galaxies for which the method extrapolates the photometric redshift, since their features lie outside the bounding box of their nearest neighbours. In these cases, the errors on photometric redshifts are larger and these galaxies are flagged appropriately.
Moreover, we analyse the impact of the galaxy density (in the feature space) on the redshift determination. In fact, there are regions in the 5-dimensional space that are more populated than others. As expected, we find that galaxies located in crowded regions have a better redshift estimation than galaxies found in sparse regions.
Since galaxies in PS1 may have incomplete photometry, our method is prepared to deal with the case of missing features. We evaluate the effect that a missing feature may produce on the results using an ablation test (artficially removing existing features). As expected, the scatter increases in these cases. This effect is especially important when the r Kron magnitude is missing, whereas a missing colour has a smaller impact.
Furthermore, we test the use of PS1 Kron colours instead of aperture colours, finding no significant difference in the results.
Although Kron colours have no physical meaning, they are easier to obtain and it is worth stressing that they can be safely used for computing photometric redshifts with our method.
Finally, we compare our results with those presented in Beck et al. (2016) for the SDSS DR12. We find that SDSS data perform slightly better than PS1 features, especially for faint galaxies (r > 20). We suggest that these differences could be caused by two main factors: different available filters, with SDSS offering a bluer band, and a stronger correlation between the SDSS magnitudes and the spectroscopic redshift. However, it is worth noticing that the overall performance of our method is fully in agreement, within the uncertainties, with the SDSS results.
A version of the code and training set is available for download at the following webpage: www.testaddress.com.