I use '-wrapdateline' option in ogr2ogr to split the polygons at the anti meridian. But in ogr2ogr Python port bWrapDateline is hardcoded to False.
Is there a way to pass 'wrapdateline' option to gdal.VectorTranslate() to split a single polygon into a multipolygon? I'm using GDAL 2.1.0
Answer
GDAL >= 2.2 supports GeoJSON RFC7946 standart. Geometries across the antimeridian are split automatically. See autotests
Sample code:
from osgeo import gdal
import json
def read_file(filename):
vsifile = gdal.VSIFOpenL(filename,'r')
gdal.VSIFSeekL(vsifile, 0, 2)
vsileng = gdal.VSIFTellL(vsifile)
gdal.VSIFSeekL(vsifile, 0, 0)
return gdal.VSIFReadL(1, vsileng, vsifile)
gdal.UseExceptions() #Fail when can't open!
gdal.VectorTranslate('/vsimem/out.json', """{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "geometry": { "type": "Polygon", "coordinates":[[[-188.56933593749997,61.58549218152362],[-173.9794921875,61.58549218152362],[-173.9794921875,66.00015035652659],[-188.56933593749997,66.00015035652659],[-188.56933593749997,61.58549218152362]]] } },
]
}""", format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])
got = read_file('/vsimem/out.json')
print json.loads(got)
Output:
{u'type': u'FeatureCollection', u'bbox': [171.4306641, 61.5854922, -173.9794922, 66.0001504], u'features': [{u'geometry': {u'type': u'MultiPolygon', u'coordinates': [[[[-180.0, 66.0001504], [-180.0, 61.5854922], [-173.9794922, 61.5854922], [-173.9794922, 66.0001504], [-180.0, 66.0001504]]], [[[180.0, 61.5854922], [180.0, 66.0001504], [171.4306641, 66.0001504], [171.4306641, 61.5854922], [180.0, 61.5854922]]]]}, u'type': u'Feature', u'properties': {}, u'bbox': [171.4306641, 61.5854922, -173.9794922, 66.0001504]}]}
No comments:
Post a Comment