[GIS] Selecting by attributes and making new shapefile using Python ogr only

Question

I have a problem making my selected features into a new shapefile with the same projection. I am new to python, any ideas?

Below is my sample script.

from osgeo import ogr
import os, sys

driver = ogr.GetDriverByName ("ESRI Shapefile")

fnVec= raw_input("Enter Filename:")
ofn=raw_input("Output Name:")
ctype=raw_input("Crime Type:")
cyear=int(raw_input("Year of incidence:"))

datasource = driver.Open (fnVec,0 )

inLayer = datasource.GetLayer()

#Execute queries(SQL)
CrimeSQL ="SELECT * FROM Crimes WHERE Crime_Type='%s' AND Year='%s'" %(ctype,cyear)
Crime =datasource.ExecuteSQL(CrimeSQL)

Crime_feature = Crime.GetNextFeature()
while Crime_feature:
    #print Crime_feature.GetField('Count')
    Crime_feature = Crime.GetNextFeature()

---

My problem now is to make the selected points into a new shapefile. I made new scripts, but I cannot view the selected points. Here is my sample code:

from osgeo import ogr, osr
import os, sys

# Get the input Layer
inDriver = ogr.GetDriverByName("ESRI Shapefile")
fnVec= raw_input("Enter Filename:")
ofn=raw_input("Output Name:")
ctype=raw_input("Crime Type:")
cyear=int(raw_input("Year of incidence:"))
inDataSource = inDriver.Open(fnVec, 0)
inLayer = inDataSource.GetLayer()
spatialRef = inLayer.GetSpatialRef()

# input SpatialReference
feature = inLayer.GetNextFeature()
geom = feature.GetGeometryRef()
inSpatialRef = geom.GetSpatialReference()

#Execute SQL
CrimeSQL = "SELECT * FROM Crimes WHERE Crime_Type='%s' AND Year='%s'" %(ctype,cyear)
Filter = inDataSource.ExecuteSQL(CrimeSQL)

# Create the output LayerS
outShapefile = os.path.join( os.path.split( fnVec )[0], ofn )
outDriver = ogr.GetDriverByName("ESRI Shapefile")

# Remove output shapefile if it already exists
if os.path.exists(ofn):
    outDriver.DeleteDataSource(ofn)

# Create the output shapefile
outDataSource = outDriver.CreateDataSource(outShapefile)
out_lyr_name = os.path.splitext( os.path.split( outShapefile )[1] )[0]
outLayer = outDataSource.CreateLayer( out_lyr_name, geom_type=ogr.wkbPoint )

# Add input Layer Fields to the output Layer if it is the one we want
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    fieldName = fieldDefn.GetName()
    outLayer.CreateField(fieldDefn)

# Get the output Layer's Feature Definition
outLayerDefn = outLayer.GetLayerDefn()


# Add features to the ouput Layer
for inFeature in inLayer:
    # Create output Feature
    outFeature = ogr.Feature(outLayerDefn)

    # Add field values from input Layer
    for i in range(0, outLayerDefn.GetFieldCount()):
        fieldDefn = outLayerDefn.GetFieldDefn(i)
        fieldName = fieldDefn.GetName()
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(),inFeature.GetField(i))
        outLayer.CreateFeature(outFeature)
print "Process Completed!"     
# Close DataSources
inDataSource.Destroy()
outDataSource.Destroy()

Answer 1

I think that with your 'CrimeSQL' request you have all desired features selected. You don't need the 'GetNextFeature' method. To test my approach, I created your fields 'Year' and 'Crime_Type' in my point shapefile of next image. There are only two features with attributes ('2011', 'Robbery') at attributes table (delimited by red rectangle).

Your code was modified by me for adapting it to the conditions of my system: Python Console of QGIS. The SQL request has 'Crime_Type = Robbery' and 'Year = 2011'.

from osgeo import ogr
import os, sys

driver = ogr.GetDriverByName ("ESRI Shapefile")

fnVec = "/home/zeito/pyqgis_data/random_points.shp"
ofn = "/home/zeito/pyqgis_data/random_points_out.shp"
ctype = "Robbery"
cyear = int("2011")

datasource = driver.Open (fnVec,0 )

inLayer = datasource.GetLayer()

#Execute queries(SQL)
CrimeSQL = "SELECT * FROM random_points WHERE Crime_Type='%s' AND Year='%s'" %(ctype,cyear)
Crime = datasource.ExecuteSQL(CrimeSQL)

multipoint = ogr.Geometry(ogr.wkbMultiPoint)

for feature in inLayer:
    geom = feature.GetGeometryRef()
    xy = geom.GetPoint()
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(xy[0], xy[1])
    multipoint.AddGeometry(point)

print multipoint.ExportToWkt()

After running the code, I used the geometry printed at the Python Console with the 'QuickWKT' plugin of QGIS to visualize the result. Features were selected as expected with your request.