Coordinate Transformation¶

In many seismic studies, we only work on a small area of few degrees, converting lat/lon to x/y coordinates is much easier for data analysis. Here are examples to use pyproj for cartographic projections.

The Mercator projection (merc) is commonly used for plotting maps.

For simple projection purpose, we will mainly use local projections, such as Stereographic projection (stere), Oblique Stereographic Alternative projection (sterea), and Gnomonic projection (gnom).

An overview of different projection views can also be found in cartopy

In [1]:

from pyproj import Proj
lat0 = 36.20
lon0 = -138.25
R = 6371.0 * 1e3 #km

from pyproj import Proj lat0 = 36.20 lon0 = -138.25 R = 6371.0 * 1e3 #km

In [2]:

proj = Proj(f"+proj=sterea +lat_0={lat0} +lon_0={lon0} +R={R} +units=km")
x,y = proj(longitude=lon0+1, latitude=lat0+1)
print(f"x = {x:.2f}, y = {y:.2f}")
lat, lon = proj(longitude=x, latitude=y, inverse=True)
print(f"lon = {lon:.2f}, lat = {lat:.2f}")

proj = Proj(f"+proj=sterea +lat_0={lat0} +lon_0={lon0} +R={R} +units=km") x,y = proj(longitude=lon0+1, latitude=lat0+1) print(f"x = {x:.2f}, y = {y:.2f}") lat, lon = proj(longitude=x, latitude=y, inverse=True) print(f"lon = {lon:.2f}, lat = {lat:.2f}")

x = 88.58, y = 111.66
lon = 37.20, lat = -137.25

In [3]:

proj = Proj(f"+proj=stere +lat_0={lat0} +lon_0={lon0} +R={R} +units=km")
x,y = proj(longitude=lon0+1, latitude=lat0+1)
print(f"x = {x:.2f}, y = {y:.2f}")
lat, lon = proj(longitude=x, latitude=y, inverse=True)
print(f"lon = {lon:.2f}, lat = {lat:.2f}")

proj = Proj(f"+proj=stere +lat_0={lat0} +lon_0={lon0} +R={R} +units=km") x,y = proj(longitude=lon0+1, latitude=lat0+1) print(f"x = {x:.2f}, y = {y:.2f}") lat, lon = proj(longitude=x, latitude=y, inverse=True) print(f"lon = {lon:.2f}, lat = {lat:.2f}")

x = 88.58, y = 111.66
lon = 37.20, lat = -137.25

In [4]:

proj = Proj(f"+proj=gnom +lat_0={lat0} +lon_0={lon0} +R={R} +units=km")
x,y = proj(longitude=lon0+1, latitude=lat0+1)
print(f"x = {x:.2f}, y = {y:.2f}")
lat, lon = proj(longitude=x, latitude=y, inverse=True)
print(f"lon = {lon:.2f}, lat = {lat:.2f}")

proj = Proj(f"+proj=gnom +lat_0={lat0} +lon_0={lon0} +R={R} +units=km") x,y = proj(longitude=lon0+1, latitude=lat0+1) print(f"x = {x:.2f}, y = {y:.2f}") lat, lon = proj(longitude=x, latitude=y, inverse=True) print(f"lon = {lon:.2f}, lat = {lat:.2f}")

x = 88.59, y = 111.67
lon = 37.20, lat = -137.25

In [5]:

proj = Proj(f"+proj=merc +lon_0={lon0} +R={R} +units=km")
x,y = proj(longitude=lon0+1, latitude=lat0+1)
print(f"x = {x:.2f}, y = {y:.2f}")
lat, lon = proj(longitude=x, latitude=y, inverse=True)
print(f"lon = {lon:.2f}, lat = {lat:.2f}")

proj = Proj(f"+proj=merc +lon_0={lon0} +R={R} +units=km") x,y = proj(longitude=lon0+1, latitude=lat0+1) print(f"x = {x:.2f}, y = {y:.2f}") lat, lon = proj(longitude=x, latitude=y, inverse=True) print(f"lon = {lon:.2f}, lat = {lat:.2f}")

x = 111.19, y = 4462.02
lon = 37.20, lat = -137.25