Thursday, 28 January 2016

python - gdal.VectorTranslate returns an empty file


When testing the following code:


from osgeo import gdal
import json


def read_file(filename):
f = gdal.VSIFOpenL(filename, "rb")
if f is None:

return None
content = gdal.VSIFReadL(1, 10000, f).decode('UTF-8')
gdal.VSIFCloseL(f)
return content


gdal.VectorTranslate('out.json', """{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "geometry": { "type": "Polygon", "coordinates": "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')
gdal.Unlink('/vsimem/out.json')

print json.loads(got)

from this answer, the file out.json is empty. The folder is with all permissions for writing allowed. When testing the method gdal.OpenEx, is possible to load an object, but still the method gdal.VectorTranslate returns an empty geojson when trying to pass this new loaded object as a srcDS object. Any hint on solving this issue?



Answer






  1. Your JSON is invalid - you have entered the "coordinates" key twice.




  2. You are trying to read from a different file than you are writing to. You are trying to write to out.json which will be a file on disk in the current working directory. You need to write to '/vsimem/out.json'




To pick up these issues more easily, use gdal.UseExceptions() so GDAL raises an exception instead of returning None when it can't open something.


I've also updated your read_file script to find the length automatically.



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!

gjs = '''{
"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]]] } },
]
}'''

ds = gdal.VectorTranslate('/vsimem/out.json', gjs, format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])


print(ds)

got = read_file('/vsimem/out.json')

#Works
print(got)

#Fails as the json module can't read you json (nor can https://jsonlint.com)
print(json.loads(got))


And once you have valid GeoJSON, gdal.VectorTranslate will still return an empty dataset. This is a known python GDAL "gotcha". The data doesn't get written until the dataset is closed and dereferenced.


This works:


from osgeo import gdal
gdal.UseExceptions()
srcDS = gdal.OpenEx('input.geojson')
ds = gdal.VectorTranslate('output.geojson', srcDS=srcDS, format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])

#Dereference and close dataset, then reopen.
del ds

ds = gdal.OpenEx('output.geojson')

print(srcDS.GetLayer(0).GetFeatureCount())
print(ds.GetLayer(0).GetFeatureCount())

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...