I'm dealing with shapefiles projected in:
GEOGCS["GCS_WGS_1984", DATUM["WGS_1984", SPHEROID["WGS_84",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]]
And I am experiencing technical difficulties to calculate a distance in (kilo)meters between two features. I can only get this distance in degrees. Here's my code:
def PointsOnLines(self, pointsLayer, linesLayer): """ Returns the number of points situated on lines+(line.width) """ # pointsLayer: shapefile's layer with point geometries # linesLaye: shapefile's layer with line geometries nbFound = 0 for i in range(pointsLayer.GetFeatureCount()): pointFeat = pointsLayer.GetNextFeature() pointGeo = pointFeat.GetGeometryRef() for j in range(linesLayer.GetFeatureCount()): lineFeat = linesLayer.GetNextFeature() try: lineGeo = lineFeat.GetGeometryRef() lineWidth = lineFeat.GetFieldAsDouble("WD1") # Gets the width of the line segment (meters) except: pass d = lineGeo.Distance(pointGeo) # Computes the shortest distance between the 2 geometries (degrees) # if d<lineWidth then point is on the line, nbFound+=1
Should I reproject the shapefile and calculate afterwards? Ans if so, which projection is best suitable (=>meters units)?
Note that I can't use this post because I don't even know the 2 points I'm currently working with.