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()