I'm automatically creating a point shapefile from a postgreSQL point table (having a geometry field called "geo" defined as geo geometry(Point,4258)) for a spesific flatenr (%%a) using
ogr2ogr --config GDAL_DATA "C:\OSGeo4W64\share\gdal" -f "ESRI Shapefile" C:\kmfelt\shape\utm32\2017\%%a_p.shp PG:"host=%Server_host% user=%Schema_user% dbname=%Database% password=sl_nab" -sql "SELECT * FROM kmfelt_o2_p where flatenr=%%a AND komid < 1800" -t_srs EPSG:%SRID% -geomfield geo
followed by some editing of the shapefile using python
call python %prog%\05_deletefields_p.py C:\kmfelt\shape\utm32\2017\%%a_p.shp
In general, this works fine and results in an output point shapefile.
However: Some times it turns out that the shapefile is empty, e.g the shapefile does not contain any points since the postgreSQL point table does not contain any data for the specific flatenr.
An empty shapefile is created, but the geometry is set to linestring?!?
... # Set POINT-shapefile fn = r""+sys.argv+"" # Read filename from input print('filnavn er ' +fn) #get the driver driver = osgeo.ogr.GetDriverByName('ESRI Shapefile') # will select the driver for our shp-file # open the data source dataSource = driver.Open(fn, 1) # 0 means read-only. 1 means writeable (read/write) # get layer layer = dataSource.GetLayer() geotype = layer.GetGeomType() print(geotype)
This results in (output on the screen)
filnavn er C:\kmfelt\shape\utm32\2017\2060_p.shp 2
The 2 indicates that the geometry type is wbkLineString. Why? (This is confirmed in QGIS: the shapefile is shown with a line symbol in the Layers panel and when loading a point style file an error message appear telling QGIS cannot apply style to layer with a different geometry type.)
How may I change/set/force the geometry type to be point (wkbPoint)?
So far I have be using this code the py-file in case of no points
... if layer.GetFeatureCount() == 0: km = ogr.FieldDefn('KMNR', ogr.OFTInteger) km.SetWidth(4) #4 integer width print "Legger til egenskapen KMNR som siste kolonne" layer.CreateField(km) # After adding fields, get the layer definition featureDefn = layer.GetLayerDefn() # Create a feature from scratch feature = ogr.Feature(featureDefn) # Populate the attribute table with a default (dummy) value print "Setter KMNR lik 9999" feature.SetField("KMNR", 9999) # Finally, feature must be written to the layer layer.CreateFeature(feature) # Clean up memory feature.Destroy()