Geographic Information Systems Asked by alxgrh on June 26, 2021
I’m trying to split a GeoJSON polygon which crosses the antimeridian. I use this solution as base gdal.VectorTranslate returns an empty file
My code is:
geojson_str = <some geojson string>
geojson = json.loads(gjstr)
print('Original GeoJSON:')
print(json.dumps(geojson, indent=2))
with open('in.geojson','w') as f:
json.dump(geojson, f)
gdal.UseExceptions()
g = gdal.OpenEx('in.geojson')
out = gdal.VectorTranslate('./out.geojson', g, format = 'GeoJSON', layerCreationOptions = ["-wrapdateline", '-lco','RFC7946=YES'])
del out
geojson_transformed = json.load(open('./out.geojson'))
print('Transformed GeoJSON:')
print(json.dumps(geojson_transformed, indent=2))
The issue is that gdal.VectorTranslate() has indeterministic behaviuor on different input.
If I use GeoJSON from example, it splits polygon into multipolygon correctly:
Original GeoJSON:
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
172.56933593749997,
61.58549218152362
],
[
186.0205078125,
61.58549218152362
],
[
186.0205078125,
66.00015035652659
],
[
172.56933593749997,
66.00015035652659
],
[
172.56933593749997,
61.58549218152362
]
]
]
}
}
]
}
Transformed GeoJSON:
{
"type": "FeatureCollection",
"name": "in",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "MultiPolygon",
"coordinates": [
[
[
[
180.0,
61.5854922
],
[
180.0,
66.0001504
],
[
172.5693359,
66.0001504
],
[
172.5693359,
61.5854922
],
[
180.0,
61.5854922
]
]
],
[
[
[
-180.0,
66.0001504
],
[
-180.0,
61.5854922
],
[
-173.9794922,
61.5854922
],
[
-173.9794922,
66.0001504
],
[
-180.0,
66.0001504
]
]
]
]
}
}
]
}
The example input contains shifted coordinates, where all negative longitudes incresed by 360. The output contains non-shifted coordinates though.
If I try to convert the example coordinates to be non-shifted, gdal.VectorTranslate() doesn’t split the polygon:
Original GeoJSON:
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
172.56933593749997,
61.58549218152362
],
[
-173.9794921875,
61.58549218152362
],
[
-173.9794921875,
66.00015035652659
],
[
172.56933593749997,
66.00015035652659
],
[
172.56933593749997,
61.58549218152362
]
]
]
}
}
]
}
Transformed GeoJSON:
{
"type": "FeatureCollection",
"name": "in",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
172.5693359,
61.5854922
],
[
172.5693359,
66.0001504
],
[
-173.9794922,
66.0001504
],
[
-173.9794922,
61.5854922
],
[
172.5693359,
61.5854922
]
]
]
}
}
]
}
So I can assume that dateline splitting is being performed only if coordinates are centered at longitude 180 (or -180).
Nevertheless when I try to use it on my custom Geojson with shifted coords, then gdal.VectorTranslate() fails to complete:
Original GeoJSON:
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
182.24833333333333,
61.336666666666666
],
[
179,
60.0
],
[
169.0,
54.0
],
[
158.9961111111111,
50.10194444444445
],
[
157.8,
52.25
],
[
157.5361111111111,
53.325833333333335
],
[
158.18333333333334,
53.79666666666667
],
[
156.75,
55.45
],
[
153.9897222222222,
58.01416666666667
],
[
155.9,
58.86666666666667
],
[
160.16666666666666,
60.583333333333336
],
[
162.21666666666667,
61.43333333333333
],
[
164.25,
62.083333333333336
],
[
174.75,
61.833333333333336
],
[
182.24833333333333,
61.336666666666666
]
]
]
},
"properties": {
"designator": "UHPP",
"color": "gray"
}
}
]
}
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-218-e81a49bf4254> in <module>()
33 gdal.UseExceptions()
34 g = gdal.OpenEx('in.geojson')
---> 35 out = gdal.VectorTranslate('./out.geojson', g, format = 'GeoJSON', layerCreationOptions = ["-wrapdateline", '-lco','RFC7946=YES'])
36 del out
37 geojson_transformed = json.load(open('./out.geojson'))
1 frames
/usr/lib/python3/dist-packages/osgeo/gdal.py in wrapper_GDALVectorTranslateDestName(*args)
3139 def wrapper_GDALVectorTranslateDestName(*args):
3140 """wrapper_GDALVectorTranslateDestName(char const * dest, Dataset srcDS, GDALVectorTranslateOptions options, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset"""
-> 3141 return _gdal.wrapper_GDALVectorTranslateDestName(*args)
3142 class GDALDEMProcessingOptions(_object):
3143 """Proxy of C++ GDALDEMProcessingOptions class."""
RuntimeError: Terminating translation prematurely after failed
translation of layer in (use -skipfailures to skip errors)
If I use non-shifted coordintates in GeoJSON input, gdal.VectorTranslate() doesn’t fail, but produce untransformed result:
Original GeoJSON:
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
-177.75166666666667,
61.336666666666666
],
[
179,
60.0
],
[
169.0,
54.0
],
[
158.9961111111111,
50.10194444444445
],
[
157.8,
52.25
],
[
157.5361111111111,
53.325833333333335
],
[
158.18333333333334,
53.79666666666667
],
[
156.75,
55.45
],
[
153.9897222222222,
58.01416666666667
],
[
155.9,
58.86666666666667
],
[
160.16666666666666,
60.583333333333336
],
[
162.21666666666667,
61.43333333333333
],
[
164.25,
62.083333333333336
],
[
174.75,
61.833333333333336
],
[
-177.75166666666667,
61.336666666666666
]
]
]
},
"properties": {
"designator": "UHPP",
"color": "gray"
}
}
]
}
Transformed GeoJSON:
{
"type": "FeatureCollection",
"name": "in",
"features": [
{
"type": "Feature",
"properties": {
"designator": "UHPP",
"color": "gray"
},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
-177.7516667,
61.3366667
],
[
174.75,
61.8333333
],
[
164.25,
62.0833333
],
[
162.2166667,
61.4333333
],
[
160.1666667,
60.5833333
],
[
155.9,
58.8666667
],
[
153.9897222,
58.0141667
],
[
156.75,
55.45
],
[
158.1833333,
53.7966667
],
[
157.5361111,
53.3258333
],
[
157.8,
52.25
],
[
158.9961111,
50.1019444
],
[
169.0,
54.0
],
[
179.0,
60.0
],
[
-177.7516667,
61.3366667
]
]
]
}
}
]
}
Is it possible to determine why gdal.VectorTranslate() fails on the third input?
Are there other reliable and precise solutions to split polygons around antimeridian? I tried geopandas-based solution from here GDAL VectorTranslate creates shards/fragments along anti-meridian, but it doesn’t seem to work.
Okay, I haven't found any suitable solution for splitting custom polygon by the antimeridian. I implemented it by myself using Pyproj and Shapely. Maybe someone will find it useful.
from pyproj import CRS, Transformer
from shapely.geometry import Point, LineString
def intersects_antimeridian_precise(start_point, end_point):
'''
Detects if line intersects the Antimeridian. Use custom Mercator projection for convenient result.
'''
# Initialize custom Mercator projection
mercator_am_proj = CRS.from_proj4("+proj=merc lon_0={}".format(start_point.x))
mercator_am_tf = Transformer.from_crs(mercator_am_proj.geodetic_crs, mercator_am_proj)
mercator_am_tf_inv = Transformer.from_crs( mercator_am_proj,mercator_am_proj.geodetic_crs)
# Project Antimeridian to the Mercator plane
antimeridian= [(180,89.99),(180,-89.99)]
am_sp_m = Point(mercator_am_tf.transform(antimeridian[0][0],antimeridian[0][1]))
am_ep_m = Point(mercator_am_tf.transform(antimeridian[1][0],antimeridian[1][1]))
antimeridian_merc = LineString([am_sp_m,am_ep_m])
# Project tested line to the Mercator plane
sp_m = Point(mercator_am_tf.transform(start_point.x, start_point.y))
ep_m = Point(mercator_am_tf.transform(end_point.x, end_point.y))
line_m = LineString([sp_m,ep_m])
return line_m.intersects(antimeridian_merc)
def wrap_antimeridian(polygon):
'''
Splits polygon which intersects the Antimeridian into two new polygons stitched to Antimeridian
Params:
polygon - list of (x,y) tuples
Returns:
List of two lists of (x,y) tuples or original polygon if it doesn't intersect the Antimeridian
'''
# Initialize Mercator transformer
mercator_am_proj = CRS.from_proj4("+proj=merc lon_0=179")
mercator_am_tf = Transformer.from_crs(mercator_am_proj.geodetic_crs, mercator_am_proj)
mercator_am_tf_inv = Transformer.from_crs( mercator_am_proj,mercator_am_proj.geodetic_crs)
antimeridian= [(180,89.99),(180,-89.99)]
am_sp_m = Point(mercator_am_tf.transform(antimeridian[0][0],antimeridian[0][1]))
am_ep_m = Point(mercator_am_tf.transform(antimeridian[1][0],antimeridian[1][1]))
antimeridian_merc = LineString([am_sp_m,am_ep_m])
forw_intersect_idx = None
forw_intersect_line = None
back_intersect_idx = None
back_intersect_line = None
# Direction in which forward line intersects the Antimeridian. Use 1 for West-East and -1 for East-West
forw_intersect_direction = None
for i,p in enumerate(polygon):
start_point = Point(p[0],p[1])
if i == len(polygon)-1:
end_point = Point(*polygon[0])
else:
end_point = Point(*polygon[i+1])
# Find intersection points
if intersects_antimeridian_precise(start_point, end_point):
if forw_intersect_idx is None:
forw_intersect_idx = i
forw_intersect_line = [start_point, end_point]
else:
back_intersect_idx = i
back_intersect_line = [start_point, end_point]
if forw_intersect_idx is not None and back_intersect_idx is not None:
# Convert lat,lon coordinates to Mercator cartesian plane and find forward and backward intersection points
sp_f_merc = Point(mercator_am_tf.transform(*forw_intersect_line[0].coords[0]))
ep_f_merc = Point(mercator_am_tf.transform(*forw_intersect_line[1].coords[0]))
forw_intersect_line_merc = LineString([sp_f_merc,ep_f_merc])
forw_intersect_point_merc = forw_intersect_line_merc.intersection(antimeridian_merc)
forw_intersect_point = mercator_am_tf_inv.transform(forw_intersect_point_merc.x,forw_intersect_point_merc.y)
sp_b_merc = Point(mercator_am_tf.transform(*back_intersect_line[0].coords[0]))
ep_b_merc = Point(mercator_am_tf.transform(*back_intersect_line[1].coords[0]))
back_intersect_line_merc = LineString([sp_b_merc,ep_b_merc])
back_intersect_point_merc = back_intersect_line_merc.intersection(antimeridian_merc)
back_intersect_point = mercator_am_tf_inv.transform(back_intersect_point_merc.x,back_intersect_point_merc.y)
# Determine intersection direction
if sp_f_merc.x < forw_intersect_point_merc.x:
forw_intersect_direction = 1
else:
forw_intersect_direction = -1
# Cut polygon before first intersection line and append first intersection point. Set longitude sign according to direction
poly_1 = polygon[:forw_intersect_idx+1]
poly_1.append((forw_intersect_point[0]*forw_intersect_direction, forw_intersect_point[1] ))
# Cut polygon after second intersection and append second intersection point before
poly_1.append((back_intersect_point[0]*forw_intersect_direction, back_intersect_point[1] ))
poly_1 += polygon[back_intersect_idx+1:]
# Cut polygon between intersection points, insert intersection points on its edges
poly_2 = [(forw_intersect_point[0]*forw_intersect_direction*-1, forw_intersect_point[1])] + polygon[forw_intersect_idx+1:back_intersect_idx+1] + [(back_intersect_point[0]*forw_intersect_direction*-1, back_intersect_point[1])]
return [poly_1, poly_2]
return polygon
Answered by alxgrh on June 26, 2021
I made a test with the json that fails for you. I mean this one
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
182.24833333333333,
61.336666666666666
],
[
179,
60.0
],
[
169.0,
54.0
],
[
158.9961111111111,
50.10194444444445
],
[
157.8,
52.25
],
[
157.5361111111111,
53.325833333333335
],
[
158.18333333333334,
53.79666666666667
],
[
156.75,
55.45
],
[
153.9897222222222,
58.01416666666667
],
[
155.9,
58.86666666666667
],
[
160.16666666666666,
60.583333333333336
],
[
162.21666666666667,
61.43333333333333
],
[
164.25,
62.083333333333336
],
[
174.75,
61.833333333333336
],
[
182.24833333333333,
61.336666666666666
]
]
]
},
"properties": {
"designator": "UHPP",
"color": "gray"
}
}
]
}
Test with ogr2ogr into GeoJSON and OpenJUMP JML formats:
ogr2ogr -f jml split.jml split.json -wrapdateline
ogr2ogr -f geojson split2.json split.json -wrapdateline -lco RFC7946=YES
In both cases the result is a multipolygon that is split by the dateline.
MULTIPOLYGON (((180.0 60.411493073371,179 60,169 54,158.996111111111 50.1019444444444,157.8 52.25,157.536111111111 53.3258333333333,158.183333333333 53.7966666666667,156.75 55.45,153.989722222222 58.0141666666667,155.9 58.8666666666667,160.166666666667 60.5833333333333,162.216666666667 61.4333333333333,164.25 62.0833333333333,174.75 61.8333333333333,180.0 61.4855893902349,180.0 60.411493073371)),((-177.751666666667 61.3366666666667,-180 60.411493073371,-180 61.4855893902349,-177.751666666667 61.3366666666667)))
I do not know what is the reason for the failure with your code but GDAL seems to be able to handle your input with ogr2ogr.
Answered by user30184 on June 26, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP