tools.projection

tools.projection#

This module provides spatial utility functions for projecting camera footprints onto the sky and matching coordinates to survey fields.

skysurvey.tools.projection.radecmodel_to_skysurface(radecmodel, ntrial=200000.0, frac=True)[source]#

Compute the sky area covered by points drawn from a ModelDAG model in RA/Dec space.

This function samples points from a ModelDAG model, projects them onto a unit sphere, and computes the convex hull of the projected points to estimate the sky area. The area can be returned as a fraction of the total sky (4π steradians) or in steradians.

Parameters:

  • radecmodel (object) – A ModelDAG-compatible model that generates RA/Dec points.

  • ntrial (int, optional) – Number of points to sample from the model. Default is 2e5. Tests suggest 2e5 is good at 0.01%.

  • frac (bool, optional) – If True, return the area as a fraction of the total sky (4π steradians). If False, return the area in steradians. Default is True.

Returns:

The sky area covered by the sampled points. If frac=True, the value is a fraction of the total sky. If frac=False, the value is in steradians.

Return type:

float

Notes

  • The RA/Dec points are converted to radians and projected using sin(dec) to account for spherical geometry.

  • The convex hull of the projected points is computed to estimate the sky area.

  • The total sky area is 4π steradians, which corresponds to 41253 square degrees.

Examples

>>> from modeldag import ModelDAG
>>> radecmodel = ...  # Your ModelDAG-compatible RA/Dec model
>>> area_frac = radecmodel_to_skysurface(radecmodel, ntrial=1e5, frac=True)
>>> print(f"Fraction of sky covered: {area_frac:.4f}")
skysurvey.tools.projection.project_to_radec(verts_or_polygon, ra, dec)[source]#

Project a geometry (or its vertices) to given ra, dec coordinates.

Parameters:

  • verts_or_polygon (shapely.Polygon or 2d-array) – geometry or vertices representing the camera footprint in the sky if vertices, the format is: x, y = vertices

  • ra (float or array) – poiting(s) R.A.

  • dec (float or array) – poiting(s) declination

Returns:

If input are vertices, returns a list of new verticies. If input are geometries, returns a list of new geometries.

Return type:

list

skysurvey.tools.projection.spatialjoin_radec_to_fields(radec, fields, how='inner', predicate='intersects', index_radec='index_radec', allow_dask=True, **kwargs)[source]#

Join the radec coordinates with the fields.

Parameters:

  • radec (DataFrame or 2d-array) –

    coordinates of the points.

    • DataFrame: must have the “ra” and “dec” columns.

      This will use the DataFrame’s index are data index.

    • 2d array (shape N,2): returned index will be ‘range(len(ra))’

  • fields ([geopandas.geoserie, geopandas.geodataframe or dict]) –

    Fields contains the fieldid and fields shapes. Several forms are accepted:

    • dict: {fieldid: 2d-array, fieldid: 2d-array …}

      here, the 2d-array are the field’s vertices.

    • geoserie: geopandas.GeoSeries with index as fieldid and geometry as field’s vertices.

    • geodataframe: geopandas.GeoDataFrame with the ‘fieldid’ column and geometry as field’s vertices.

Returns:

(geometry.sjoin result)

Return type:

GeoDataFrame

skysurvey.tools.projection.parse_fields(fields)[source]#

Read various formats for fields and returns it as a geodataframe.

Parameters:

fields ([geopandas.geoserie, geopandas.geodataframe or dict]) –

Fields contains the fieldid and fields shapes. Several forms are accepted:

  • dict: {fieldid: 2d-array or regions, fieldid: 2d-array or regions …}

    here, the 2d-array are the field’s vertices or a astropy/ds9 regions

  • geoserie: geopandas.GeoSeries with index as fieldid and geometry as field’s vertices.

  • geodataframe: geopandas.GeoDataFrame with the ‘fieldid’ column and geometry as field’s vertices.

Return type:

GeoDataFrame (geometry.sjoin result)

Examples

provide a dict of ds9 regions >>> fields = {450:”box(50,30, 3,4,0)”, 541:”ellipse(190,-10,1.5,1,50)”} >>> geodf = parse_fields(fields)

skysurvey.tools.projection.regions_to_shapely(region)[source]#

Converts astropy Region into a shapely geometry.

Parameters:

region (str or Regions (see astropy-regions.readthedocs.io)) –

  • If str, it is assumed to be in the dr9 ircs format, e.g.: region = "box(40.0, 50.0, 5.0, 4.0, 0.0)".

  • If Regions, region will be converted into the str format, using region = region.serialize("ds9").strip().split("\n")[-1].

The following format have been implemented:

  • box

  • circle

  • ellipse

  • polygon

Returns:

the geometry will depend on the input regions.

Return type:

Shapely’s Geometry

Raises:

NotImplementedError – if the format is not recognised.

Examples

>>> shapely_ellipse = regions_to_shapely('ellipse(54,43.4, 4, 2,-10)')
>>> shapely_rotated_rectangle = regions_to_shapely('box(-30,0.4, 4, 2,80)')
skysurvey.tools.projection.cart2sph(vec)[source]#

Converts cartesian [x,y,z] to spherical [r, theta, phi] coordinates (in degrees).

Parameters:

vec (array) – x, y, z

Returns:

[r, theta, phi]

Return type:

array

skysurvey.tools.projection.sph2cart(vec)[source]#

Converts spherical coordinates [r, theta, phi] to cartesian coordinates [x,y,z].

Parameters:

vec (array) – r, theta, phi ; angles in degrees

Returns:

[x, y, z]

Return type:

array

skysurvey.tools.projection.rot_xz(vec, theta)[source]#

Rotates cartesian vector v [x,y,z] by angle theta around axis (0,1,0).

Parameters:

  • vec (array) – x, y, z

  • theta (float) – angle in degree

Returns:

rotated x, y, z

Return type:

array

skysurvey.tools.projection.rot_xz_sph(l, b, theta)[source]#

Rotate spherical coordinate (l,b = theta, phi) by angle theta around axis (0,1,0) (calls does to rot_xz() and cart2sph()).

Parameters:

  • l (float) – spherical coordinate

  • b (float) – spherical coordinate

  • theta (float) – angle in degree

Returns:

[r, theta, phi]

Return type:

array