You make a mistake with the pyproj parameters
from pyproj import Proj, transform
inProj = Proj("+init=EPSG:3857"))
outProj = Proj("+init=EPSG:4326")
x1,y1 = (-20037507.0672, 20037509.6184)
x2,y2 = (-1467048.29156, 8625918.8737)
# with the x,y coordinates of pointx
print transform(inProj,outProj,x1,y1)
(-179.99998854118687, 85.08398750388278)
# with the x,y coordinates of point pointy
print transform(inProj,outProj,x2,y2)
(-13.178719028497799, 61.16289280460126)
Control with GDAL (Python GDAL/OGR Cookbook: Projections)
from osgeo import osr
from_srs = osr.SpatialReference()
from_srs.ImportFromEPSG(3857)
to_srs = osr.SpatialReference()
to_srs.ImportFromEPSG(4326)
transf = osr.CoordinateTransformation(from_srs,to_srs)
print transf.TransformPoint(x1,y1)[:2]
(-179.99998854118684, 85.05112976833757)
# with the x,y coordinates of point y
print transf.TransformPoint(x2,y2)[:2]
(-13.178719028497795, 61.000416666720845)
If you want to use lists of coordinates:
x = (-20037507.0672,-1467048.29156)
y = (20037509.6184, 8625918.8737)
for i,j in zip(x,y):
print transform(inProj,outProj,i, j)
(-179.99998854118687, 85.08398750388278)
(-13.178719028497799, 61.16289280460126)
This should do the trick. I haven't been able to install the pyproj module, but at the point marked #do transformation. I tested it by doing simple arithmetic. I put in the function based on your code.
Basic structure is:
- Load data from file.
- Traverse json string - essentially nested dictionaries and lists
- Iterate through the coordinates. The coordinate lists are heavily nested to allow for multipart polygons
- Write altered data to a new file
Code:
import pyproj
import json
#define projections
wgs84 = Proj('epsg:4326')
IrishGrid = Proj('epsg:29902')
#load in data
with open(r'path_to_file_you_want_to_edit.json', 'r') as f:
data = json.load(f)
#traverse data in json string
for feature in data['features']:
#print feature['geometry']['type']
#print feature['geometry']['coordinates']
#all coordinates
coords = feature['geometry']['coordinates']
#coordList is for each individual polygon
for coordList in coords:
#each point in list
for coordPair in coordList:
print coordPair
x1 = coordPair[0]
y1 = coordPair[1]
lat_grid, lon_grid = numpy.meshgrid(x1, y1)
#do transformation
coordPair[0],coordPair[1] = pyproj.transform(IrishGrid,wgs84,lat_grid, long_grid)
#write reprojected json to new file
with open('path_to_new_file.json', 'w') as f:
f.write(json.dumps(data))
Hope this helps.
Best Answer