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

Compare your code with the cookbook example https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-shapefile-and-add-data. — user30184, Feb 14 '19 at 10:17
This is what I did, I took my code from "Convert polygon to points" https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html. — ilFonta, Feb 14 '19 at 10:22
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? — user30184, Feb 14 '19 at 10:35
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. — ilFonta, Feb 14 '19 at 10:49
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 https://www.gdal.org /ogr__core_8h.html#a800236a0d460ef66e687b7b65610f12a. — user30184, Feb 14 '19 at 11:21
Hi, I didn't modified anything from the example, I'm not so skilled. I followed your suggestion, I changed point.AddPoint(Xcoord, Ycoord) in point.AddPoint_2D(Xcoord, Ycoord). Now I don't have z-coordinate, but it still doesn't works. I think the bug is here outFeature.SetGeometry(multipoint) — ilFonta, Feb 14 '19 at 18:20

score 2 · Answer 1 · answered Feb 19 '19 at 17:59

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()

Empty shapefile created with Python and GDAL

1 Answers1