docs: def query(self, coords, order=1):
          """
          Returns the map value at the specified location(s) on the sky.
  
          Args:
              coords (`astropy.coordinates.SkyCoord`): The coordinates to query.
              order (Optional[int]): Interpolation order to use. Defaults to `1`,
                  for linear interpolation.
  
          Returns:
              A float array containing the map value at every input coordinate.
              The shape of the output will be the same as the shape of the
              coordinates stored by `coords`.
          """
          out = np.full(len(coords.l.deg), np.nan, dtype='f4')
  
          for pole in self.poles:
              m = (coords.b.deg >= 0) if pole == 'ngp' else (coords.b.deg < 0)
  
              if np.any(m):
                  data, w = self._data[pole]
                  x, y = w.wcs_world2pix(coords.l.deg[m], coords.b.deg[m], 0)
                  out[m] = map_coordinates(data, [y, x], order=order, mode='nearest')
  
          return out