Set or change the geometry type of an empty point shapefile using Python with GDAL/OGR

Question

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.

Example: flatenr = 2060 does not exist in the table

An empty shapefile is created, but the geometry is set to linestring?!?

From 05_deletefields.py

...
# Set POINT-shapefile
fn = r""+sys.argv[1]+""   # 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()   

Answer 1

The solution turned out to be quite simple. I had deeper look into the options of ogr2ogr and discovered the

-lco NAME=VALUE: Layer creation option (format specific)

Since I'm creating av shapefile, I had a look at GDAL's ESRI Shapefile format page. In the section "Layer Creation Options" I found

SHPT=type: Override the type of shapefile created. Can be one of NULL for a simple .dbf file with no .shp file, POINT, ARC, POLYGON ...

So adding -lco SHPT=POINT in my ogr2ogr call, did the trick:

ogr2ogr --config GDAL_DATA "C:\OSGeo4W64\share\gdal" -f "ESRI Shapefile" C:\kmfelt\shape\utm32\2017\%%a_p.shp -lco SHPT=POINT 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

An empty POINT shapefile is created.