1

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 pcjericks.github.io/py-gdalogr-cookbook/…. – user30184 Feb 14 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 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 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 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 at 11:21
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() 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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