0

I am converting the geographic coordinates used on old historic maps into WGS84. Through educated guesswork I can establish all the necessary parameters to rebuild the grid system and generate a wkt file for the Coordinate Reference System.

Is it possible to used this customised CRS file with pyproj and how would I go about referencing the file and converting the data?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Tymon
  • 11
  • 2
  • Not sure I understand the problem. You have your WKT in a file, so to re-use it you read the file in (using standard python open/read) and then feed the text to the PyProj transformation methods. Do you know how to use PyProj with any standard CRS? – Spacedman Apr 02 '22 at 12:34

2 Answers2

1

Given some WKT in a text file in my case in /tmp/mycrs.wkt:

PROJCRS["OSGB 1936 / British National Grid",BASEGEOGCRS["OSGB 1936",DATUM["OSGB 1936",ELLIPSOID["Airy 1830",6377563.396,299.3249646,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4277]],CONVERSION["British National Grid",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",49,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996012717,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",400000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",-100000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["unknown"],AREA["UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E"],BBOX[49.75,-9.2,61.14,2.88]],ID["EPSG",27700]]

I read it in using standard python I/O and create a projection object:

mycrs = open("/tmp/mycrs.wkt","r").read()
p = pyproj.Proj(mycrs)

and then I can use that to convert from map coords to anything, say lat-long epsg:4326:

In [33]: pyproj.transform(p, "epsg:4326", 500000,600000)
Out[33]: (55.283987628697844, -0.42711506505796637)
Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • Mind updating your example to use CRS & Transformer? https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 – snowman2 Apr 02 '22 at 15:48
1

Thanks that was just the guidance I needed. I upgraded the code:

#get custom Coordinate Reference System.
mycrs = open(r"GSGS_3868_Ed_4.prj", "r").read()
inproj = pyproj.crs.CRS(mycrs)

#set output CRS
outproj = 2326


#create transformation
proj = pyproj.Transformer.from_crs(inproj, outproj, always_xy=True)


calculate new location


x2,y2 = proj.transform(val_x,val_y)
Tymon
  • 11
  • 2