0

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:

Distances from evaluation

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.

  • 4
    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. – FSimardGIS Oct 16 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 at 20:53
  • 1
    Have a look at epsg.io for New York, it lists a series of projections for the state. – FSimardGIS Oct 16 at 20:55
  • 1
    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. – FSimardGIS Oct 16 at 21:01
  • @FSimardGIS you can convert your comments into an answer. – MarcM Oct 17 at 9:58
2

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.

| improve this answer | |
  • 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. – snowman2 Oct 18 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! – FSimardGIS Oct 18 at 4:17
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')]
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.