2

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.

enter image description here

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()
  • 1
    Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format. You had an ArcPy tag on your question but make no mention of using that software in your question and so I have removed it. – PolyGeo May 29 '16 at 6:30
  • Please do not include "Thanks" in your posts here. These Q&As are not personal communications between users but instead are wiki-style posts and I'll wager that you have never seen "Thanks" at the bottom of a Wikipedia page. – PolyGeo May 29 '16 at 6:48
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).

enter image description here

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.

enter image description here

  • Thanks, 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: – Chris Odan Jun 1 '16 at 0:36

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.