This is a follow up to a question from Chris Jackson (Imperial) on LinkedIn.
Here are the coordinates we have, UTM zone 16, WGS84 datum:
x = 812136.38
y = 10151410.68
Those look like they are in feet. Let's convert them:
x *= 0.3048
y *= 0.3048
x, y
(247539.168624, 3094149.9752640002)
We will use the pyproj
project.
import pyproj as pp
Set up the projections using EPSG codes (we'll define the target projection too while we're at it).
utm15_wgs84 = pp.Proj(init='epsg:32615')
utm16_wgs84 = pp.Proj(init='epsg:32616')
utm15_nad27 = pp.Proj(init='epsg:26715')
utm16_nad27 = pp.Proj(init='epsg:26716')
Let's just check where we are with a quick map:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
% matplotlib inline
--------------------------------------------------------------------------- ImportError Traceback (most recent call last) <ipython-input-5-a31ec74bf0f1> in <module>() ----> 1 from mpl_toolkits.basemap import Basemap 2 import matplotlib.pyplot as plt 3 get_ipython().magic('matplotlib inline') ImportError: No module named 'mpl_toolkits.basemap'
m = Basemap(width=6000000,height=6000000,projection='lcc',
resolution='c',
lat_1=45,
lat_2=55,
lat_0=50,
lon_0=-100.)
# Set up the map.
m.drawcoastlines()
m.drawmapboundary(fill_color='#aaccee')
m.fillcontinents(color='#aaddaa',lake_color='#aaccee')
lon, lat = utm16_wgs84(x, y, inverse=True) # To convert to lon, lat
map_x, map_y = m(lon, lat)
m.scatter(map_x, map_y, 12, marker='o', color='red')
plt.show()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-6-0b773f659872> in <module>() ----> 1 m = Basemap(width=6000000,height=6000000,projection='lcc', 2 resolution='c', 3 lat_1=45, 4 lat_2=55, 5 lat_0=50, NameError: name 'Basemap' is not defined
Looks reasonable, let's do the transformation...
new_coords = pp.transform(utm16_nad27, # From
utm15_wgs84, # To
x, y) # Coordinates
new_coords # These are in metres!
([668289.7894208364, 668389.5137875683, 668489.2381542905], [498.66323992501503, 598.3968759964243, 698.1308414315228])
Sanity check with another map:
m = Basemap(width=6000000,height=6000000,projection='lcc',
resolution='c',
lat_1=45,
lat_2=55,
lat_0=50,
lon_0=-100.)
# Set up the map.
m.drawcoastlines()
m.drawmapboundary(fill_color='#aaccee')
m.fillcontinents(color='#aaddaa',lake_color='#aaccee')
new_x, new_y = new_coords
new_lon, new_lat = utm15_nad27(new_x, new_y, inverse=True) # To convert to lon, lat
new_map_x, new_map_y = m(new_lon, new_lat)
m.scatter(new_map_x, new_map_y, 12, marker='o', color='red')
plt.show()
Sweet!
© 2015 Agile Geoscience — CC-BY — Have fun!