I’ve got a series of points and polygons represented in lon/lat pairs. Each set of points that I’ll be working with lie within 3km of each other, but I have different clusters from different points in the US. I’m trying to convert them to ‘map XY’ coordinates in meters, so that I can do simply geometric/distance calculations: point is within polygon, distance between points, distance from point to polygon, etc.
I’ve attempted something with pyproj using a transformer: pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857") but I’m clearly doing something wrong. When using the converted XY pairs, the calculated distance is the same for (-90., 30.) to (-90.5, 30.) compared with (-90., 50.) to (-90.5, 50.), which can’t be correct as the second pair of points is much closer to the north pole.
The projection you are converting to, “EPSG:3857” is pseudo-Mercator. It is the flattened map used by Google Maps etc. It is not great at all for doing distance computations at large distances, as indeed the distance between longitudes 90 and 90.5 does not depend on latitude in this projection, as two red lines that correspond to your two lines show.
The only good part about Mercator for distance computations, and the reason to convert from 4326, is it preserves shapes of the geometries. So around any given point, the “stretch” is equal in horizontal and vertical directions. Thus you can compute a local stretch coefficient for a given location (set of points close enough to each other) for converting 3857 distances to real distances, and for a small enough location it gives you reasonable results. But the stretch coefficient would vary between distant locations, like Ontario and Louisiana in your examples.
You can also get better precision by using different local projections for each location / set of points.