Why are the 500 meters not 500 meters? And how do I add 500 meters in x/y direction to a lat/lon point?

Question

I wanted to create a square cell on a map with size 500m times 500m. Therefore, I wrote the following code:

from pyproj import Transformer

bl_lat, bl_lon = 40.64423086189233, -74.0307935018417
trans_coord_to_meter = Transformer.from_crs('epsg:4326', 'epsg:3857')
bl_t = trans_coord_to_meter.transform(bl_lat, bl_lon)

print(trans_coord_to_meter.transform(bl_t[0], bl_t[1], direction='INVERSE'))  # blue
# 40.64423086189233, -74.0307935018417 (lat, lon)
print(trans_coord_to_meter.transform(bl_t[0]+500, bl_t[1], direction='INVERSE'))  # red, + x_axis
# 40.64423086189233, -74.02630192542111
print(trans_coord_to_meter.transform(bl_t[0], bl_t[1]-500, direction='INVERSE'))  # green, - y_axis
# 40.640822707278815, -74.0307935018417
print(trans_coord_to_meter.transform(bl_t[0]+500, bl_t[1]-500, direction='INVERSE'))  # yellow, + x and y axis

I thought that this adds 500m in both dimensions (maybe I misunderstood the function). But, when I check the result, I observed the following:

So, how do I really get points an actual distance of 500m in between? (for x and y)

I want to keep using Python, but I'm flexible with the library.

Answer 1

Beware of EPSG:3857 / Mercator projection for measurements... It distorts distances a lot, especially when far from the equator. 500 Mercator meters is only about 500 meters times cos(latitude) in reality. (380 meters in your case at 40°).

Instead, you should use a projection better suited for your area of interest. Some websites, like epsg.io, can display a list of coordinate systems for a specified region. For the United States, there is also the option of checking which State Plane coordinate system correspond to your area of interest. NGS NCAT tool can automatically determine the State Plane coordinate system for a specified Lat/Lon location and convert it.

According to NCAT, your area falls in NY L-3104 New York Long Island (epsg 32118 for the NAD83 crs). That crs could be worth a try.

Answer 2

Building on the answer by @FSimardGIS:

In pyproj 3, you will be able to query the PROJ database for CRS in a specific location: https://pyproj4.github.io/pyproj/latest/api/database.html#pyproj-database-query-crs-info

You can try the dev version of pyproj 3: https://pypi.org/project/pyproj/3.0.dev3/

>>> from pyproj.database import query_crs_info, query_utm_crs_info
>>> from pyproj.aoi import AreaOfInterest
>>> bl_lat, bl_lon = 40.64423086189233, -74.0307935018417
>>> aoi = AreaOfInterest(bl_lon, bl_lat, bl_lon, bl_lat)
>>> crs_list = query_crs_info(auth_name="EPSG", pj_types="PROJECTED_CRS", area_of_interest=aoi)
>>> len(crs_list)
46
>>> pprint([(crs.name, crs.code) for crs in crs_list if "NAD83(2011)" in crs.name])
[('NAD83(2011) / UTM zone 18N', '6347'),
 ('NAD83(2011) / Conus Albers', '6350'),
 ('NAD83(2011) / New Jersey', '6526'),
 ('NAD83(2011) / New Jersey (ftUS)', '6527'),
 ('NAD83(2011) / New York Long Island', '6538'),
 ('NAD83(2011) / New York Long Island (ftUS)', '6539')]
>>> pprint([(crs.name, crs.code) for crs in crs_list if "NAD83 " in crs.name])
[('NAD83 / New York Long Island (ftUS)', '2263'),
 ('NAD83 / UTM zone 18N', '26918'),
 ('NAD83 / New Jersey', '32111'),
 ('NAD83 / New York Long Island', '32118'),
 ('NAD83 / Statistics Canada Lambert', '3347'),
 ('NAD83 / New Jersey (ftUS)', '3424'),
 ('NAD83 / Canada Atlas Lambert', '3978'),
 ('NAD83 / BLM 18N (ftUS)', '4438'),
 ('NAD83 / Conus Albers', '5070')]
>>> utm_crs_list = query_utm_crs_info(area_of_interest=aoi)
>>> pprint([(crs.name, crs.code) for crs in utm_crs_list])
[('NAD27 / UTM zone 18N', '26718'),
 ('NAD83 / UTM zone 18N', '26918'),
 ('WGS 72 / UTM zone 18N', '32218'),
 ('WGS 72BE / UTM zone 18N', '32418'),
 ('WGS 84 / UTM zone 18N', '32618'),
 ('NAD83(NSRS2007) / UTM zone 18N', '3725'),
 ('NAD83(HARN) / UTM zone 18N', '3748'),
 ('NAD83(2011) / UTM zone 18N', '6347')]