Empty shapefile created with Python and GDAL

Question

I need to create a point shapefile using the coordinates specified in the variable I called "my_points_dataset".

'my_points_dataset' variable has this structure:

1615311.055743243079633 4820229.073873873800039 0,1615321.055743243079633 4820229.073873873800039 0,1615331.055743243079633 4820229.073873873800039 0,1615341.055743243079633 4820229.073873873800039 0,1615351.055743243079633 4820229.073873873800039 0, etc`

The class of this variable is 'osgeo.ogr.Geometry'

Each point in 'my_points_dataset' is defined as x, y and z coordinates (z could be negligeable because is always 0). x, y and z data are separated by a space, each point is separated by comma.

I wrote this Python code, but at the end the result is an empty shapefile. Do yo know why?

import ogr, gdal
shpDriver = ogr.GetDriverByName('ESRI Shapefile')
 if os.path.exists('/path/mypoints.shp'):
 shpDriver.DeleteDataSource('/path/points_shapefile.shp')
 outDataSource = shpDriver.CreateDataSource('/path/points_shapefile.shp')
 outLayer = outDataSource.CreateLayer('/path/points_shapefile.shp', geom_type=ogr.wkbMultiPoint)
 featureDefn = outLayer.GetLayerDefn()
 outFeature = ogr.Feature(featureDefn)
 outFeature.SetGeometry(my_points_dataset)
 outLayer.CreateFeature(outFeature)
 outFeature = None

• user30184: Compare your code with the cookbook example pcjericks.github.io/py-gdalogr-cookbook/….
• ilFonta: This is what I did, I took my code from "Convert polygon to points" pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html.
• user30184: It seem that you have skipped many steps before doing outFeature.SetGeometry(my_points_dataset). Have you tried to use the code of the tutorial exactly?
• ilFonta: Yes I did. I copied/pasted the example from pcjerks in spyder and I used a general shapefile with only one polygon as test. The firsta part of the example works and the point net is created and saved in my_points-dataset. The conversion to shapefile doesn't work.
• user30184: It seems that you have modified the example by adding z coordinate, is that right? I guess that wkbMultiPoint is not correct geometry type for your x, y, z data. Drop z coordinates and try again. I am not sure what type matches your data but perhaps it is wkbMultiPoint25D gdal.org /ogr__core_8h.html#a800236a0d460ef66e687b7b65610f12a.

Answer 1

I found the solution in gene's answer in Writing a shapefile using OGR, from Shapely Geometry - No feature added Error

 # Write point coordinates to Shapefile
shpDriver = ogr.GetDriverByName('ESRI Shapefile')
# Set the crs
latlong = osr.SpatialReference()
latlong.ImportFromEPSG( 3003 )
if os.path.exists('/path/mypoints.shp'):
 shpDriver.DeleteDataSource('/path/mypoints.shp')
outDataSource = shpDriver.CreateDataSource('/path/mypoints.shp')
outLayer = outDataSource.CreateLayer('', srs=latlong, geom_type=ogr.wkbMultiPoint)
outLayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
outDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(outDefn)
outFeature.SetField('id', 1)
outFeature.SetGeometry(my_points_dataset)
outLayer.CreateFeature(outFeature)
# Remove temporary files
outDataSource.Destroy()