2

I am working on a bunch of shapefile data, and I need to process it and after some kind of processing, I need to feed the data to an indexer (in this case ElasticSearch). The problem is I just discovered that there are some anomalies in the data. I found this in one of my polygons.

[(593.0587471450001, -1.8046608309999783),
 (593.0587471160001, -1.8046608379999611),
 (593.0587470810001, -1.804660872999932),
....]

Now correct me if I am wrong, but aren't longitudes suposed to be in the range of -180 to 180 and latitudes in the range of -90 to 90? So the numbers above just don't make sense. Is there any function in shapely, GeoPandas or any other library that can validate/check whether a coordinate is valid or not? Or should I do that manually?

UPDATE:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Also I don't know if this is relevant, but I only found this one polygon with very big numbers in it. The other polygons in the same file are fine (in the right range).

2
  • @user2856 That's the first time I have ever heard of projected coordinates. I am so new in this. Thanks. I updated my question with the information from the .crs attribute. – catris25 May 28 '20 at 3:38
  • It might help to think about how these coordinates were generated. If they were produced by some automated process, it's possible they're just not properly wrapped and you might be able to correct it with a simple script, e.g., by subtracting 360. – Rob Skelly Jun 28 '20 at 20:40
1

use crs attribute, for example

x = gpd.read_file("your_shapefile.shp")
coord = x.crs
print(coord)

there should be an EPSG value, and that is your coordinate system.
according to this right here what you need to do about your over 180 degrees coordinate value is to subtract/add it by 360 until it is in the range of -180 and 180 degrees.

to verify if your data is correct, you might want another shapefile, that you know that the shapefile is in good projection quality, as your orientation and overlay/subplot them.

#read trusted good quality shapefile for orientation
verif = gpd.read_file("orientation_shapefile.shp")

#plot the layers
base = x.plot()
verif.plot(ax=base)
plot.show()
0

Standing on the shoulders of the Q and several commentaries, this might do the trick. Note that it could return a boolean result, but since it has a side-effect of wrapping the longitude, it returns the modified (lon,lat) tuple, if it is in the bounds of the WGS84.


def valid_lonlat(lon: float, lat: float) -> Optional[Tuple[float, float]]:
    """
    This validates a lat and lon point can be located
    in the bounds of the WGS84 CRS, after wrapping the
    longitude value within [-180, 180)

    :param lon: a longitude value
    :param lat: a latitude value
    :return: (lon, lat) if valid, None otherwise
    """
    # Put the longitude in the range of [0,360):
    lon %= 360
    # Put the longitude in the range of [-180,180):
    if lon >= 180:
        lon -= 360
    lon_lat_point = shapely.geometry.Point(lon, lat)
    lon_lat_bounds = shapely.geometry.Polygon.from_bounds(
        xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0
    )
    # return lon_lat_bounds.intersects(lon_lat_point)
    # would not provide any corrected values

    if lon_lat_bounds.intersects(lon_lat_point):
        return lon, lat

# it would be used like:

lon_lat = valid_lonlat(lon, lat)
if lon_lat is None:
    raise ValueError("(lon, lat) is not in WGS84 bounds")
else:
    lon, lat = lon_lat

For a possible OOP approach (imports omitted):


LON_LAT_BOUNDS = shapely.geometry.Polygon.from_bounds(
    xmin=-180.0, ymin=-90.0, xmax=180.0, ymax=90.0
)


@dataclass_json
@dataclass(order=True)
class LonLat:
    """Coordinates for longitude and latitude (assumed to be in degrees)"""

    #: longitude (in degrees)
    lon: float
    #: latitude (in degrees)
    lat: float
    #: units (degrees)
    units: str = "degrees"

    lon_lat_bounds = LON_LAT_BOUNDS

    def __post_init__(self):
        # Put the longitude in the range of [0,360):
        self.lon %= 360
        # Put the longitude in the range of [-180,180):
        if self.lon >= 180:
            self.lon -= 360

    @property
    def point(self) -> shapely.geometry.Point:
        return shapely.geometry.Point(self.lon, self.lat)

    def is_valid(self) -> bool:
        """
        This validates a lat and lon point can be located
        in the bounds of the WGS84 CRS, after wrapping the
        longitude value within [-180, 180)

        :return: True if the (lon,lat) is in WGS84 bounds
        """
        return LON_LAT_BOUNDS.intersects(self.point)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.