1

I have a .shp file with a large number of coordinates and I'm not sure how to transform to GPS coordinates. Ideally I'd like to use Python to convert , but have not seen this formatting before:

[(4738040.670931757, 949345.8379265075),
(4738040.04068242, 949233.7326115519),
(4737188.151903, 948466.4097770005),]

I know the location is Oregon USA.

6
  • If I am right, in qgis you can export the shp to gpx, is that the output format you are looking for?
    – Gery
    Nov 13, 2021 at 0:45
  • Gpx? Do you mean GPS?
    – Derek_P
    Nov 13, 2021 at 1:36
  • I know qgis provides gpx as an export option, not sure if gps (ie. file format extension, *.gps) is also provided by qgis as an export format option.
    – Gery
    Nov 13, 2021 at 2:02
  • Have a look at: pcjericks.github.io/py-gdalogr-cookbook/projection.html#
    – Ben W
    Nov 13, 2021 at 2:35
  • I'll have to check that out - thanks!
    – Derek_P
    Nov 13, 2021 at 2:41

1 Answer 1

2

By GPS coordinates, I guess you mean WGS84. A geographic, geocentric Coordinate Reference System, using angular latitude and longitude coordinates, with degrees as units. It looks like your coordinates are in a projected coordinate system using (I guess) feet as units.

Most spatial Python libraries have some methods for transforming between Coordinate Reference Systems.

To simply transform your list of coordinate pairs you can use the Shapely and PyProj libraries. I guessed at a source CRS for your data and used NAD83(HARN)/Oregon South epsg:2914, and converted to WGS84, epsg:4326 with the following code snippet:

from shapely.geometry import Point
from shapely.ops import transform
import pyproj

source = pyproj.CRS('EPSG:2914')
target = pyproj.CRS('EPSG:4326')


project = pyproj.Transformer.from_crs(source, target, always_xy=True).transform
    
coord_list = [(4738040.670931757, 949345.8379265075),
              (4738040.04068242, 949233.7326115519),
              (4737188.151903, 948466.4097770005)]

transformed_coord_list = []

for coord_pair in coord_list:
    point = Point(coord_pair[0], coord_pair[1])
    transformed_coord = transform(project, point)
    transformed_coord_list.append((transformed_coord.y, transformed_coord.x))

print('Transformed coords (lat, long order): ', transformed_coord_list)

The output from the print statement is:

Transformed coords (lat, long order):  [(44.26929990402178, -121.19940821027379), (44.26899241512979, -121.19940704224689), (44.26686829750889, -121.20263417672915)]

On the other hand, you could use Geopandas to transform your entire shapefile to WGS84.

import geopandas as gpd

in_path = 'C:\\Path\\to\\input.shp'

out_path = 'C:\\Path\\to\\output_WGS84.shp'

in_file = gpd.read_file(in_path)

out_file_WGS84 = in_file.to_crs({'init': 'epsg:4326'})

out_file_WGS84.to_file(out_path)

I created a test point shapefile from your coordinates in epsg:2914 and converted to epsg:4326 (WGS84) using the code above. You can see the result below loaded into QGIS over Open Street Map and labelled with geometry values (in lat, long order), which plots the points in Redmond, Oregon.

enter image description here

You can find Shapely docs on transformations here.

Geopandas docs here.

Examples of reprojecting vector data with Geopandas here.

You could also use the ogr and osr modules which are included in the GDAL Python bindings to do the same tasks. There are some useful recipes to be found here.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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