Python GeoJSON – Convert Geometry to New Coordinate System and WKB Format

geojsonpythonshapelywell-known-binary

I have a GeoJSON like such:

{
    "type": "FeatureCollection",
    "features": [{
        "type": "Feature",
        "properties": {
            "uuid": "68F15C01-20CD-4D77-954B-9483BA1D4D33",
            "name": "ABC",
            "len": "529",
            "infotype": "Restricted",
            "h1": null,
            "join": null
        },
        "geometry": {
            "type": "LineString",
            "coordinates": [
                [-360909.60698310119, 7600968.922204642, 0.0],
                [-361357.344715965, 7600811.951385159, 0.0],
                [-361805.08159795138, 7600654.939420643, 0.0]
            ]
        }
    }]
}

I need to convert the geometry to EPSG:27700, as I've just realised I have no idea what the coordinate system above is supposed to be but I do know the placement of the features should be in the UK.

Once I've successfully done that, I'd then like the different feature geometries in WKB format using Python.
I've tried the below code but I'm getting the error:

ValueError: could not convert string to float: '-360909.60698310117, 7600968.922204642, 0.0, -361357.344715965, 7600811.951385159, 0.0, -361805.08159795136, 7600654.939420643, 0.0'

# Python program to read
# json file


import json
import shapely
from shapely import wkb, wkt
from shapely.geometry import Point

# Opening JSON file
f = open('metadata.json')


data = json.load(f)


for i in range(len(data['features'])):
    geom = str(data['features'][i]['geometry']['coordinates'])
    newstr = geom.replace("[", "")
    newstr2 = newstr.replace("]","")
    geom3 = float(newstr2)
    g = wkb.dumps(geom3, hex=True, srid=4326)

    print(g)

Best Answer

The coordinate system of GeoJSON may be EPSG:3857. You can use the following script

import json
import shapely
from shapely import wkb, wkt
from shapely.geometry import LineString

# Opening JSON file
f = open('metadata.json')
data = json.load(f)

for feat in data["features"]:
    geom = LineString(feat['geometry']['coordinates'])
    
    g = wkb.dumps(geom, hex=True, srid=3857)
    
    print(g)

# OUT:
# 01020000A0110F00000300000063FA8C6D360716C19E66053BD2FE5C410000000000000000D038FD60350E16C1947EE3FCAAFE5C410000000000000000D2698E53341516C1C3771FBC83FE5C410000000000000000

Related Solutions

[GIS] Transform GeoJSON geometry to WKB

Use ST_GeomFromGeoJSON.

SELECT ST_GeomFromGeoJSON('{"type":"MultiPolygon","coordinates":[[[[-125045.48351212002,4577567.588141698],[-124816.19981552364,4577552.93014355],[-124765.99472517562,4577419.175847012],[-124842.47121534991,4577392.905406596],[-125045.48351212002,4577567.588141698]]]]}');

Returns:

010600000001000000010300000001000000050000006B3477BC5787FEC0141DA4E54776514134C371320379FEC0D078873B447651411DF264EADF75FEC0D31341CB22765141541B198AA77AFEC0822EF2391C7651416B3477BC5787FEC0141DA4E547765141

To set the SRID (3857: Google Mercator) and transform to EPSG 4326:

ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON('{"type":"MultiPolygon","coordinates":[[ ... ]]]}'), 3857), 4326)

Note that without setting the SRID, PostGIS will "assign" the SRID of 0. Your question did not indicate a projection, but it looks like Google Mercator (EPSG 3857). Your comments indicate you want 4326, so once the SRID has been correctly set, you can transform to that.

This returns:

0106000020E6100000010000000103000000010000000500000074A6F23D0CF9F1BF443A04EE50FE424051E7FD7F9CF0F1BFA74F72874DFE424040E59997C3EEF1BFE64B6F7F2EFE4240A669C2F693F1F1BF8356276728FE424074A6F23D0CF9F1BF443A04EE50FE4240

Which in WKT is:

MULTIPOLYGON(((-1.12330269049053 37.9868447800295,-1.12124300000001 37.986741,-1.12079200000001 37.985794,-1.121479 37.985608,-1.12330269049053 37.9868447800295)))

Python OGR to Shapely – How to Convert Geometries with M-Values to Those Without

Here's a potential solution following suggestions from jbalk's and user30184's comments. It uses OGR's SetMeasured function to drop the M-values and relies on WKB (instead of WKT) to do the actual conversion between OGR and Shapely.

def from_ogr_to_shapely(ogr_geom):
    # Creating a copy of the input OGR geometry. This is done in order to 
    # ensure that when we drop the M-values, we are only doing so in a 
    # local copy of the geometry, not in the actual original geometry. 
    ogr_geom_copy = ogr.CreateGeometryFromWkb(ogr_geom.ExportToIsoWkb())
    
    # Dropping the M-values
    ogr_geom_copy.SetMeasured(False)

    # Generating a new shapely geometry
    # Note: I needed to typecast the results from
    # `ogr_geom_copy.ExportToIsoWkb()` to bytes since, in newer versions 
    # of GDAL/OGR and shapely, the outputs of that method were incompatible 
    # with the expected inputs for the `shapely.wkb.loads` method, 
    # leading to an error.
    shapely_geom = shapely.wkb.loads(bytes(ogr_geom_copy.ExportToIsoWkb()))

    return shapely_geom

As mentioned in the commented section in the code above, we need to create a local copy of the input OGR geometry. This is because the function used to drop the M-values (the SetMeasured function) actually modifies the original input geometry. So in order to ensure that our original geometry remains intact, we need to create a local copy. Then we can modify the local copy to our heart's content without fear of changes creeping upstream to the original OGR geometry.

Here's the function being applied to the inputs in the question:

print(from_ogr_to_shapely(pt_2D).wkt)
# POINT (0 0) # right!

print(from_ogr_to_shapely(line1_2D).wkt)
# LINESTRING (0 0, 1 1) # right!

print(from_ogr_to_shapely(mline_2D).wkt)
# MULTILINESTRING ((0 0, 1 1), (3 1, 4 1)) # right!

print(from_ogr_to_shapely(pt_3D).wkt)
# POINT Z (0 0 0) # right!

print(from_ogr_to_shapely(line1_3D).wkt)
# LINESTRING Z (0 0 0, 1 1 1) # right!

print(from_ogr_to_shapely(mline_3D).wkt)
# MULTILINESTRING Z ((0 0 0, 1 1 1), (3 1 3, 4 1 6)) # right!

Edit from the future

As pointed out by @ermaure in their answer, it seems like more recent versions of GDAL/OGR and shapely have changed things a bit and now the outputs of the ogr_geom.ExportToIsoWkb method are incompatible with the expected inputs for the shapely.wkb.loads method, causing the shapely.wkb.loads(ogr_geom_copy.ExportToIsoWkb()) statement to yield an error. This can be fixed by forcing the result of the ExportToIsoWkb method to be typecast to bytes. Explicitly, the command needs to change to shapely.wkb.loads(bytes(ogr_geom_copy.ExportToIsoWkb())). I have edited the code above to reflect that change.

Special thanks to @ermaure. Please be sure to upvote their answer as well!

Best Answer

Related Solutions

[GIS] Transform GeoJSON geometry to WKB

Python OGR to Shapely – How to Convert Geometries with M-Values to Those Without

Edit from the future

Related Question