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

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.

• Beware of Mercator projection for measurements... It does not preserve distances, especially when far from the equator. 500 Mercator meters is only about 500 meters times cos(latitude) (380 meters in your case at 40°). Instead you should use a projection better suited for your area of interest. Oct 16, 2020 at 20:49
• @FSimardGIS Thank you! Where can I find an overview of these projections? Or maybe be the one suitable for New York City?
– So S
Oct 16, 2020 at 20:53
• Have a look at epsg.io for New York, it lists a series of projections for the state. Oct 16, 2020 at 20:55
• For instance, according to NGS NCAT, that area falls in NY L-3104 New York Long Island (epsg 32118 for the NAD83 crs). That crs could be worth a try. Oct 16, 2020 at 21:01

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.

• This was a very helpful answer. As a side note: epsg,io can be out of date and is not the official source for EPSG data: github.com/pyproj4/pyproj/issues/711#issuecomment-705807124. epsg.org is the official source. Oct 18, 2020 at 0:06
• @snowman2 Thanks for the comment. Indeed, epsg.org is the official registry, somehow I have this reflex of searching on epsg.io, but I must remember to double-check the info :-) But now that the registry's got a new website recently it's great for searching and retrieving the projection and transformation info, and you can do a map search of crs as well. I will certainly use it more often in the future! Oct 18, 2020 at 4:17

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) / 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 York Long Island', '32118'),
('NAD83 / New Jersey (ftUS)', '3424'),