Skip to content

srs

Functions:

  • centeredLAEA –

    Load a Lambert-Azimuthal-Equal_Area spatial reference system (SRS) centered

  • loadSRS –

    Load a spatial reference system (SRS) from various sources.

  • tileIndexAt –

    Get the "slippery tile" index at the given zoom, around the

  • xyTransform –

    Transform xy points between coordinate systems.

centeredLAEA 

centeredLAEA(
    lon=None, lat=None, name="unnamed_m", geom=None
)

Load a Lambert-Azimuthal-Equal_Area spatial reference system (SRS) centered on a given set of latitude and longitude coordinates. Alternatively, a geom can be passed to center the LAEA on.

Parameters:

  • lon

    (float, default: None ) –

    The longitude of the projection's center. Required if no geom is given.

  • lat

    (float, default: None ) –

    The latitude of the projection's center. Required if no geom is given.

  • geom

    –

    The region shape to center the LAEA in. If given, lat and lon must not be given, instead they will be defined automatically as the coordinates of the region centroid.

Returns:

  • SpatialReference –
Source code in geokit/core/srs.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def centeredLAEA(lon=None, lat=None, name="unnamed_m", geom=None):
    """
    Load a Lambert-Azimuthal-Equal_Area spatial reference system (SRS) centered
    on a given set of latitude and longitude coordinates. Alternatively, a geom
    can be passed to center the LAEA on.

    Parameters
    ----------
    lon : float
        The longitude of the projection's center. Required if no geom is given.

    lat : float
        The latitude of the projection's center. Required if no geom is given.

    geom: osgeo.ogr.Geometry
        The region shape to center the LAEA in. If given, lat and lon must not
        be given, instead they will be defined automatically as the coordinates
        of the region centroid.

    Returns
    -------
    osr.SpatialReference
    """
    if geom is None:
        assert isinstance(lat, numbers.Number) and isinstance(lon, numbers.Number), (
            "If geom is not passed, lat and lon must be given as float values."
        )
    else:
        assert isinstance(geom, ogr.Geometry), "geom must be given as osgeo.ogr.Geometry class object if not None."
        assert lat is None and lon is None, "If geom is given, lat and lon must not be given."

    # check if lat/lon can be used or if it must be extracted from geom first
    if not geom is None:
        # transform to EPSG:4326 in case it not already is lat/lon projection
        geom = GEOM.transform(geom, toSRS=4326)
        # extract lat/lon centroid coordinates to center LAEA upon
        lon = geom.Centroid().GetX()
        lat = geom.Centroid().GetY()

    srs = osr.SpatialReference()
    srs.ImportFromWkt(
        'PROJCS["{}",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",{}],PARAMETER["longitude_of_center",{}],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'.format(
            name, float(lat), float(lon)
        )
    )

    if gdal.__version__ >= "3.0.0":
        srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    # assert that the srs is valid (may be invalid if e.g. wrong integer codes were passed)
    assert srs.Validate() == 0, f"Created srs is invalid."

    return srs

loadSRS 

loadSRS(
    source: srs_input,
    geom: Geometry | None = None,
    **kwargs,
) -> SpatialReference

Load a spatial reference system (SRS) from various sources.

Parameters:

  • source

    (str | int | SpatialReference | None) –

    The SRS to load. Can be: * osr.SpatialReference object * EPSG integer ID * Standardized SRS string like 'EPSG:4326' * Common SRS name in SRSCOMMON (e.g. 'europe_laea') * WKT string * 'laea' — creates a Lambert Azimuthal Equal Area projection (optionally centered on a geometry)

  • geom

    (Geometry, default: None ) –

    Geometry to center the projection on if 'laea' is used.

  • **kwargs

    –

    Passed directly to SRSCOMMON.laea(**kwargs), allowing customization such as lon, lat, or name.

Returns:

  • SpatialReference –
Source code in geokit/core/srs.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def loadSRS(source: srs_input, geom: ogr.Geometry | None = None, **kwargs) -> osr.SpatialReference:
    """
    Load a spatial reference system (SRS) from various sources.

    Parameters
    ----------
    source : str | int | osr.SpatialReference | None
        The SRS to load. Can be:
          * osr.SpatialReference object
          * EPSG integer ID
          * Standardized SRS string like 'EPSG:4326'
          * Common SRS name in SRSCOMMON (e.g. 'europe_laea')
          * WKT string
          * 'laea' — creates a Lambert Azimuthal Equal Area projection
            (optionally centered on a geometry)

    geom : osgeo.ogr.Geometry, optional
        Geometry to center the projection on if 'laea' is used.

    **kwargs :
        Passed directly to `SRSCOMMON.laea(**kwargs)`, allowing
        customization such as `lon`, `lat`, or `name`.

    Returns
    -------
    osr.SpatialReference
    """
    # Do initial check of source
    if isinstance(source, osr.SpatialReference):
        return source
    elif source is None:
        return None

    # Create an empty SRS
    srs = osr.SpatialReference()

    # Check if source is a string
    if isinstance(source, str):
        src_upper = source.strip().upper()

        # Handle special LAEA case
        if src_upper == "LAEA":
            # merge geom into kwargs if provided
            if geom is not None and "geom" not in kwargs:
                kwargs["geom"] = geom
            return SRSCOMMON.laea(**kwargs)

        elif hasattr(SRSCOMMON, source):
            # assume a name for one of the common SRS's was given
            srs = SRSCOMMON[source]
        else:
            try:
                # try handling as a standardized epsg or esri etc. code
                srs = osr.SpatialReference()
                _val = srs.SetFromUserInput(source)
                assert _val == 0
            except:
                srs.ImportFromWkt(source)  # assume a Wkt string was input
    elif isinstance(source, int):
        srs.ImportFromEPSG(source)
    else:
        raise GeoKitSRSError("Unknown srs source type: ", type(source))

    if gdal.__version__ >= "3.0.0":
        srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    # assert that the srs is valid (may be invalid if e.g. wrong integer codes were passed)
    if srs.Validate() != 0:
        raise RuntimeError(f"Created srs is invalid.")

    return srs

tileIndexAt 

tileIndexAt(x, y, zoom, srs)

Get the "slippery tile" index at the given zoom, around the coordinates ('x', 'y') within the specified 'srs'.

Source code in geokit/core/srs.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
def tileIndexAt(x, y, zoom, srs):
    """Get the "slippery tile" index at the given zoom, around the
    coordinates ('x', 'y') within the specified 'srs'.
    """
    srs = loadSRS(srs)

    iterable_input = isinstance(x, Iterable) or isinstance(y, Iterable)

    if not srs.IsSame(EPSG4326):
        pt = xyTransform(x, y, fromSRS=srs, toSRS=EPSG4326, outputFormat="xy")

        if iterable_input:
            x, y = pt.x, pt.y
        else:
            x, y = pt.x[0], pt.y[0]

    if iterable_input:
        x = np.array(x)
        y = np.array(y)

    xi, yi = smopy.deg2num(y, x, zoom)

    return Tile(xi, yi, zoom)

xyTransform 

xyTransform(
    *args,
    toSRS: srs_input,
    fromSRS: srs_input,
    outputFormat: Literal["raw", "xy", "xyz"] = "raw",
) -> (
    list[tuple] | TransformedPointsXY | TransformedPointsXYZ
)

Transform xy points between coordinate systems.

Returns:

  • list of tuples, or namedtuple –
    • See the point for the 'outputFormat' argument
Source code in geokit/core/srs.py 
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def xyTransform(
    *args, toSRS: srs_input, fromSRS: srs_input, outputFormat: Literal["raw", "xy", "xyz"] = "raw"
) -> list[tuple] | TransformedPointsXY | TransformedPointsXYZ:
    """Transform xy points between coordinate systems.

    Parameters
    ----------
        xy : A single, or an iterable of (x,y) tuples
            The coordinates to transform

        toSRS : Anything acceptable by geokit.srs.loadSRS
            The srs of the output points

        fromSRS : Anything acceptable by geokit.srs.loadSRS
            The srs of the input points

        outputFormat : str
            Determine return value format
            * if 'raw', the raw output from osr.TransformPoints is given
            * if 'xy', or 'xyz' the points are given as named tuples

    Returns
    -------
    list of tuples, or namedtuple
      * See the point for the 'outputFormat' argument
    """
    # load srs's
    fromSRS = loadSRS(fromSRS)
    toSRS = loadSRS(toSRS)

    # make a transformer
    trx = osr.CoordinateTransformation(fromSRS, toSRS)

    # Do transformation
    if len(args) == 0:
        raise GeoKitSRSError("no positional inputs given")
    elif len(args) == 1:
        xy = args[0]
        if isinstance(xy, tuple):
            x, y = xy
            out = [
                trx.TransformPoint(x, y),
            ]
        else:
            out = trx.TransformPoints(xy)
    elif len(args) == 2:
        x = np.array(args[0])
        y = np.array(args[1])
        xy = np.column_stack([x, y])

        out = trx.TransformPoints(xy)

    else:
        raise GeoKitSRSError("Too many positional inputs")
    # Done!
    if outputFormat == "raw":
        return out
    elif outputFormat == "xy":
        x = np.array([o[0] for o in out])
        y = np.array([o[1] for o in out])

        return TransformedPointsXY(x, y)

    elif outputFormat == "xyz":
        x = out[:, 0]
        y = out[:, 1]
        z = out[:, 2]

        return TransformedPointsXYZ(x, y, z)