[GIS] How to validate longitude and latitude coordinates in Python

geopandaspythonshapely

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).

Best Answer

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()

Related Solutions

GDAL Grid – How to Create a Grid of 10degx10deg Polygons for a Shapefile Using GDAL

This will do it (in R):

xmin <- -180
ymin <- -90
side <- 10 # degrees of the side of each cell
listP <- list() # creates an empty list (which is similar to an array)
lands <- vector() # creates an empty array
k <- 1
nx <- 36 # number of cols
ny <- 18 # number of rows
for (i in 1:nx) { # for each col
  x1 <- xmin + (i-1)*side
  x2 <- xmin + i*side
  print(paste(i,"/",nx,sep="")) # just in case it takes too long
  flush.console() # same above
  for (j in 1:ny) { # for each row
    y1 <- ymin + (j-1)*side
    y2 <- ymin + j*side
    p = Polygon(cbind(c(x1,x1,x2,x2,x1),c(y1,y2,y2,y1,y1)))
    pP = Polygons(list(p), as.character(k))
    listP[[length(listP)+1]] <- pP # append to the array
    lands[length(lands)+1] <- 0 # append to the array
    k <- k+1
  }
}
Grid <- SpatialPolygons(listP, 1:(k-1))
d <- data.frame(LAND=lands)
spdf <- SpatialPolygonsDataFrame(Grid, data=d)

Looks like the code in Java you showed should be something like:

form osgeo import gdal
from osgeo import ogr
from osgeo import osr
import os
import random

# Specify that you want to work with a Shapefile
DriverName = "ESRI Shapefile"
# Not sure what this line does
driver = ogr.GetDriverByName(DriverName)

# Name of shapefile to create
FileName = 'gride.shp'

# simple check to see if that file already exists in the system, if so delete it
if os.path.exists(FileName):
  driver.DeleteDataSource(FileName)

# Now create an Empty shapefile Data source
shapeData = driver.CreateDataSource(FileName)

fieldName = 'LAND'

lng_min = -180
lng_max = 180
step = 10

lat_min = -90
lng_max = 90

poly = org.Geometry(ogr.wkbPolygon)
for lng_var in range(lng_min, lng_max, step):
    for lat_var in range(lat_min, lat_max, step):
            ring = org.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(lng_var, lat_var)
            ring.AddPoint(lng_var + step, lat_var)
            ring.AddPoint(lng_var + step, lat_var + step)
            ring.AddPoint(lng_var, lat_var + step)
            ring.AddPoint(lng_var, lat_var)
            poly.AddGeometry(ring)

Python – Extracting Longitude and Latitude from Shapefile

Your soa_shape_map data is in a projected coordinate reference system (CRS) that uses metres for horizontal coordinate units (# TM65_Irish_Grid). You either need to project it to the same geographic CRS as your lon, lats (I don't know what geographic CRS you're using) or project your lon, lats to the same projected CRS as your shapefile.

E.g. (assuming WGS84 as the geographic CRS)

import geopandas as gpd

soa_shape_map_path = "path/to/SA2011_Esri_Shapefile_0/SA2011.shp"
soa_shape_map = gpd.read_file(soa_shape_map_path)

soa_shape_map_geo = soa_shape_map.to_crs(epsg=4326)  # EPSG 4326 = WGS84 = https://epsg.io/4326
print(soa_shape_map.iloc[0].geometry.centroid.y, soa_shape_map_geo.iloc[0].geometry.centroid.y)

377176.1212192121 54.62955604723917

Also, to avoid your explicit nested loops, you could have a look at geopandas spatial joins, e.g:

import geopandas as gpd

points_path = r"points.shp"
points_shape_map = gpd.read_file(points_path)
points_shape_map.to_crs(epsg=4326, inplace=True)

soa_shape_map_path = r"path/to/SA2011_Esri_Shapefile_0/SA2011.shp"
soa_shape_map = gpd.read_file(soa_shape_map_path)
soa_shape_map.to_crs(epsg=4326, inplace=True)

points_polys = gpd.sjoin(points_shape_map, soa_shape_map)  # Inner join (the default), only contains intersecting points (with the attributes of intersecting poly)
# points_polys = gpd.sjoin(points_shape_map, soa_shape_map, how="inner")  # Explicitly request an inner join
points_polys = gpd.sjoin(points_shape_map, soa_shape_map, how="left")  # Outer left join, contains all points and intersecting poly attributes (leave poly attribute None or NaN where there was no poly to intersect the point)

Best Answer

Related Solutions

GDAL Grid – How to Create a Grid of 10degx10deg Polygons for a Shapefile Using GDAL

Python – Extracting Longitude and Latitude from Shapefile

Related Question