Python GDAL/OGR: WGS84 & meters units

Question

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.

Answer 1

To get distances in metres, you have to use a projection for your data. As long as it is not continent-wide, you can choose a suitable UTM projection. It does not matter if the two points are not in the same zone, but in neighbouring ones. The error is not too big if you choose one UTM zone for both (compared to calculations on a sphere).

If the data is continent-wide, you better choose a lambert projection (laea or lcc).