sketchingpy.geo

  1"""Utilities for geospatial operations.
  2
  3Utilities for geospatial operations, inspired by processing-geopoint under the BSD license at
  4https://github.com/SchmidtDSE/processing-geopoint/tree/main.
  5
  6License:
  7    BSD
  8"""
  9import itertools
 10import math
 11import typing
 12
 13import sketchingpy.shape_struct
 14
 15BASE_SCALE = 0.0001
 16BASE_LATITUDE = 0
 17BASE_LONGITUDE = 0
 18BASE_X = 0
 19BASE_Y = 0
 20
 21RADIUS_MAJOR = 6378137.0
 22RADIUS_MINOR = 6356752.3142
 23
 24TRANSFORM_MAYBE = typing.Optional['GeoTransformation']
 25
 26
 27class PixelOffset:
 28    """Structure to place a map projection into pixel space."""
 29
 30    def __init__(self, x: float, y: float):
 31        """Create PixelOffset.
 32
 33        Make data structure indicating where the map should be centered inside of a sketch in terms
 34        of pixels.
 35
 36        Args:
 37            x: Horizontal coordinate at which to center in pixels.
 38            y: Vertical coordinate at which to center in pixels.
 39        """
 40        self._x = x
 41        self._y = y
 42
 43    def get_x(self) -> float:
 44        """Get the x center position.
 45
 46        Returns:
 47            Horizontal coordinate at which to center in pixels.
 48        """
 49        return self._x
 50
 51    def get_y(self) -> float:
 52        """Get the y center position.
 53
 54        Returns:
 55            Vertical coordinate at which to center in pixels.
 56        """
 57        return self._y
 58
 59
 60class GeoPoint:
 61    """Utility to convert latitude and longitud to x and y."""
 62
 63    def __init__(self, longitude: float, latitude: float):
 64        """Create a new geographic point.
 65
 66        Create utility to convert latitude and longitude points to pixel coordinates (x and y) using
 67        a Web Mercador projection.
 68
 69        Args:
 70            longitude: The longitude of the point in degrees.
 71            latitude: The latitude of the point in degrees.
 72        """
 73        self._longitude = longitude
 74        self._latitude = latitude
 75        self._x = self._longitude_to_x(self._longitude)
 76        self._y = self._latitude_to_y(self._latitude)
 77
 78    def get_x(self, transform: TRANSFORM_MAYBE = None) -> float:
 79        """Get the horizontal pixel coordinate for this geographic point.
 80
 81        Args:
 82            transform: Map view (transformation) to use in finding this coordinate.
 83
 84        Returns:
 85            Converted x position in pixels.
 86        """
 87        if transform is None:
 88            return self._x * BASE_SCALE
 89
 90        geo_offset_x = transform.get_geo_offset().get_x()
 91        pixel_offset_x = transform.get_pixel_offset().get_x()
 92        scale = transform.get_scale()
 93        return (self._x * BASE_SCALE - geo_offset_x) * scale + pixel_offset_x
 94
 95    def get_y(self, transform: TRANSFORM_MAYBE = None) -> float:
 96        """Get the vertical pixel coordinate for this geographic point.
 97
 98        Args:
 99            transform: Map view (transformation) to use in finding this coordinate.
100
101        Returns:
102            Converted y position in pixels.
103        """
104        if transform is None:
105            return self._y * BASE_SCALE
106
107        geo_offset_y = transform.get_geo_offset().get_y()
108        pixel_offset_y = transform.get_pixel_offset().get_y()
109        scale = transform.get_scale()
110        return (self._y * BASE_SCALE - geo_offset_y) * scale + pixel_offset_y
111
112    def get_longitude(self) -> float:
113        """Get the longitude of this point.
114
115        Returns:
116            The longitude of this point in degrees.
117        """
118        return self._longitude
119
120    def get_latitude(self) -> float:
121        """Get the latitude of this point.
122
123        Returns:
124            The latitude of this point in degrees.
125        """
126        return self._latitude
127
128    def _longitude_to_x(self, longitude: float) -> float:
129        return math.radians(longitude) * RADIUS_MAJOR
130
131    def _latitude_to_y(self, latitude: float) -> float:
132        if latitude > 89.9:
133            latitude = 89.9
134        elif latitude < -89.9:
135            latitude = -89.9
136
137        return -1.0 * math.log(math.tan(
138            math.pi / 4.0 + math.radians(latitude) / 2.0
139        )) * RADIUS_MAJOR
140
141
142class GeoPolygon:
143    """Collection of GeoPoints."""
144
145    def __init__(self, points: typing.List[GeoPoint]):
146        """Create a new polygon make up of GeoPoints.
147
148        Args:
149            points: The points in this closed polygon.
150        """
151        self._points = points
152
153    def get_points(self) -> typing.List[GeoPoint]:
154        """Get the points in this polygon.
155
156        Returns:
157            List of points in this geographic polygon.
158        """
159        return self._points
160
161    def to_shape(self, transform: TRANSFORM_MAYBE = None) -> sketchingpy.shape_struct.Shape:
162        """Convert this GeoPolygon to a closed shape in pixel-space.
163
164        Args:
165            transform: The transformation or "map view" to use in finding pixel coordinates.
166
167        Returns:
168            Closed shape after projection to pixels.
169        """
170        def make_tuple(target: GeoPoint) -> typing.Tuple[float, float]:
171            return (target.get_x(transform=transform), target.get_y(transform=transform))
172
173        flat_tuples = [make_tuple(point) for point in self._points]
174
175        shape = sketchingpy.shape_struct.Shape(flat_tuples[0][0], flat_tuples[0][1])
176        for point in flat_tuples[1:]:
177            shape.add_line_to(point[0], point[1])
178        shape.close()
179
180        return shape
181
182
183class GeoTransformation:
184    """Description of a geographic transfomration or "map view" to use in projections."""
185
186    def __init__(self, geo_offset: GeoPoint, pixel_offset: PixelOffset, scale: float):
187        """Create a new map view record.
188
189        Args:
190            geo_offset: The center of the map view in terms of longitude and latitude.
191            pixel_offset: Where to place the map view in the sketch in terms of pixels.
192            scale: The map zoom level.
193        """
194        self._geo_offset = geo_offset
195        self._pixel_offset = pixel_offset
196        self._scale = scale
197
198    def get_geo_offset(self) -> GeoPoint:
199        """Get the map center in terms of longitude and latitude.
200
201        Returns:
202            The center of the map view in terms of longitude and latitude.
203        """
204        return self._geo_offset
205
206    def get_pixel_offset(self) -> PixelOffset:
207        """Get the center of this map within the sketch's pixel coordinates.
208
209        Returns:
210            Where to place the map view in the sketch in terms of pixels.
211        """
212        return self._pixel_offset
213
214    def get_scale(self) -> float:
215        """Get the map scale factor.
216
217        Returns:
218            The map zoom level.
219        """
220        return self._scale
221
222
223class GeoPolygonBuilder:
224    """Builder to create Shapes by way of GeoPolygons and GeoTransformations."""
225
226    def __init__(self, longitude: float, latitude: float,
227        transform_getter: typing.Callable[[], TRANSFORM_MAYBE]):
228        """Start a builder with a single point.
229
230        Args:
231            longitude: The longitude of the starting point in degrees.
232            latitude: The latitude of the starting point in degrees.
233            transform_getter: Function to call to get the transform to use when converting to shape.
234        """
235        self._points = [GeoPoint(longitude, latitude)]
236        self._transform_getter = transform_getter
237
238    def add_coordinate(self, longitude: float, latitude: float):
239        """Add a geographic point to this polygon.
240
241        Draw a line from the last point to a new position where this new point is defined
242        geographically in degrees.
243
244        Args:
245            longitude: The longitude of the next point in degrees.
246            latitude: The latitude of the next point in degrees.
247        """
248        self._points.append(GeoPoint(longitude, latitude))
249
250    def to_shape(self) -> sketchingpy.shape_struct.Shape:
251        """Convert this GeoPolygon to a closed shape in pixel-space.
252
253        Convert this GeoPolygon to a closed shape in pixel-space using the current map view /
254        geospatial transformation configuration.
255
256        Returns:
257            Closed shape after projection to pixels.
258        """
259        polygon = GeoPolygon(self._points)
260        transform = self._transform_getter()
261        return polygon.to_shape(transform)
262
263
264def parse_geojson(source: typing.Dict) -> typing.List[typing.List[typing.Tuple[float, float]]]:
265    """Utility to parse GeoJSON into a series of GeoPolygons.
266
267    Utility to parse GeoJSON into a series of GeoPolygons which currently only supports MultiPolygon
268    and Polygon.
269
270    Args:
271        source: The loaded GeoJSON source to parse.
272
273    Returns:
274        Points of polygons found withing the GeoJSON source.
275    """
276    all_shapes = []
277
278    for feature in source['features']:
279        geometry = feature['geometry']
280
281        feature_shapes: typing.Iterable = []
282        if geometry['type'] == 'MultiPolygon':
283            feature_shapes = itertools.chain(*geometry['coordinates'])
284        elif geometry['type'] == 'Polygon':
285            feature_shapes = [geometry['coordinates'][0]]
286        else:
287            raise RuntimeError('Only support for MultiPolygon and Polygon.')
288
289        for shape in feature_shapes:
290            points = [(float(x[0]), float(x[1])) for x in shape]
291            all_shapes.append(points)
292
293    return all_shapes