2

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
Improve this question
6
  • Compare your code with the cookbook example pcjericks.github.io/py-gdalogr-cookbook/….
    – user30184
    Feb 14, 2019 at 10:17
  • This is what I did, I took my code from "Convert polygon to points" pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html.
    – ilFonta
    Feb 14, 2019 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, 2019 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, 2019 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 gdal.org /ogr__core_8h.html#a800236a0d460ef66e687b7b65610f12a.
    – user30184
    Feb 14, 2019 at 11:21

1 Answer 1

Reset to default
2

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()
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.