1

With a Python script, I want to split a vector layer by attributes. For each unique value of a define field, I want to create a new layer with a name that contains this attribute in camelCase and with features concerned. I divide my logic in three steps :

  1. Iterate throught my source layer to find every unique value of my field of interest (city_code).
  2. Create each necessary new layer, with the good name (city_name in camelCase).
  3. Insert values in this respectives layers.

OGR give me good results for steps 1 and 2, but I can't achieve step 3. I can't iterate a second time in my features and insert values in corresponding layers I've created in step 2. Maybe my logic is wrong.

Documentation says "You must call ResetReading if you want to start iterating over the layer again.", perfect for my purpose but I'm not able to apply this.

What the best way to split a vector layer by attributes with OGR ?

Data (parcel.shp)

id | city_name | city_code
--------------------------
1  | New-York  | 001
2  | New-York  | 001
3  | New-York  | 001
4  | Toronto   | 002
5  | Toronto   | 002
6  | Toronto   | 002
7  | Seattle   | 003
8  | Seattle   | 003
9  | Seattle   | 003

Script

from osgeo import ogr
from collections import OrderedDict

# Input
file = r"parcel"
inputFolder = r"C:/Users/me/myFolder/"
outputFolder = r"C:/Users/me/myFolder/shapes/"

# Layer
extension = r".shp"
inputFile = inputFolder + file + extension
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(inputFile)
layer = dataSource.GetLayer()
spatialReference = layer.GetSpatialRef()

# Step 1

# List by attributes
listeCityCode = []
for feat in layer:
    cityCode = feat.GetField("city_code")
    listeCityCode.append(cityCode)    

listeCityCodeUnique = set(listeCityCode)
listeCityCodeUniqueOrdered = list(OrderedDict.fromkeys(listeCityCodeUnique))

# Selection by city_code
index = 0
for i in listeCityCodeUniqueOrdered:    
    sql = f'SELECT * FROM {file} WHERE city_code LIKE \'{listeCityCodeUniqueOrdered[index]}\''
    index += 1
    layerQuery = dataSource.ExecuteSQL(sql)

    # Get fields    
    for feature in layerQuery :
        id = feature.GetField("id")
        city_name = feature.GetField("city_name")
        city_code = feature.GetField("city_code")
        geom = feature.GetGeometryRef()

    # City name in camelCase
    c0 = str(city_name)
    c1 = c0.replace("-"," ")
    c2 = c1.title()
    city_name_camelCase = c2.replace(" ","")

    # Step 2

    # Create a layer by city name       
    newDataSource = driver.CreateDataSource(r"{}parcel_{}_{}.shp".format(outputFolder, city_code, city_name_camelCase))
    newLayer = newDataSource.CreateLayer("parcel_{}_{}".format(city_code, city_name_camelCase), spatialReference, ogr.wkbPolygon) 

    # Create fields
    newLayer.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
    field_city_name = ogr.FieldDefn("city_name", ogr.OFTString)
    field_city_name.SetWidth(50)
    newLayer.CreateField(field_city_name)
    field_city_code = ogr.FieldDefn("city_code", ogr.OFTString)
    field_city_code.SetWidth(10)
    newLayer.CreateField(field_city_code) 

    # Step 3

    for feature in layerQuery:         

        newFeature = ogr.Feature(newLayer.GetLayerDefn())

        newFeature.SetField("id", id)
        newFeature.SetField("city_name", city_name)
        newFeature.SetField("city_code", code_code)
        newFeature.SetGeometry(geom)           

        newLayer.CreateFeature(newFeature)

        del newDataSource, newLayer
0

Just a feedback, because I found a way to execute my problem. One of the important difficulties here is to create layers without knowing the name in advance. So I split the creation of layers and the insertion of elements into them in 2 steps. Another important point is the opening and closing of data sources. And also the read or writing mode when opening a layer.

Here is a full script, with user feedbacks :

# -*- coding: utf-8 -*-

from osgeo import ogr
from collections import OrderedDict
import glob
import os

# Input settings
myFile = r"parcel"
inputFolder = r"C:/Users/me/myFolder/"
outputFolder = r"C:/Users/me/myFolder/mySubFolder/"

# Source layer
extension = r".shp"
inputFile = inputFolder + myFile + extension
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(inputFile)
layer = dataSource.GetLayer()
spatialReference = layer.GetSpatialRef()

print(f"---- Beginning of process ----\n\nLayer analysis :\n{myFile}{extension}")

# Step 1 : List of city codes in source layer
listCityCode = []
for feat in layer:
    cityCode = feat.GetField("city_code")
    listCityCode.append(cityCode)    

listCityCodeUnique = set(listCityCode)
listCityCodeUniqueOrdered = list(OrderedDict.fromkeys(listCityCodeUnique))
numberCityCode = len(listCityCodeUniqueOrdered)

print(f"|-> {numberCityCode} city codes found\n\nWork in progress, please wait...")

# Step 2 : create a shapefile by city code

# Selection of city code in source layer
index = 0
for i in listCityCodeUniqueOrdered:    
    sql = f'SELECT * FROM {myFile} WHERE city_code LIKE \'{listCityCodeUniqueOrdered[index]}\''
    index += 1
    layerQuery = dataSource.ExecuteSQL(sql)

    # Get attributes    
    for feature in layerQuery :        
        city_name = feature.GetField("city_name")
        city_code = feature.GetField("city_code")

    # CamelCase and specials characters (especially for non-English languages)
    c0 = str(city_name)
    c1 = c0.replace("-"," ").replace("'"," ").replace("à","a").replace("â","a").replace("ä","a").replace("é","e").replace("è","e").replace("ë","e").replace("ê","e").replace("í","i").replace("ì","i").replace("ï","i").replace("î","i").replace("ô","o").replace("ö","o").replace("ó","o").replace("ò","o").replace("ú","u").replace("ù","u").replace("û","u").replace("ü","u").replace("æ","ae").replace("œ","oe").replace("ç","c").replace("ñ","n")
    c2 = c1.title().replace("À","A").replace("Â","A").replace("Ä","A").replace("É","E").replace("È","E").replace("Ê","E").replace("Ë","E").replace("Í","I").replace("Ì","I").replace("Ï","I").replace("Î","I").replace("Ô","O").replace("Ö","O").replace("Ó","O").replace("Ò","O").replace("Ú","U").replace("Ù","U").replace("Û","U").replace("Ü","U").replace("Æ","Ae").replace("Œ","Oe").replace("Ç","C").replace("Ñ","N")
    city_name_camelCase = c2.replace(" ","")

    # Create a layer by city code        
    newDataSource = driver.CreateDataSource(r"{}parcel_{}_{}.shp".format(outputFolder, city_code, city_name_camelCase))
    newLayer = newDataSource.CreateLayer("{}parcel_{}_{}".format(outputFolder, city_code, city_name_camelCase), spatialReference, ogr.wkbPolygon) 

    # Create fields
    field_id = ogr.FieldDefn("id", ogr.OFTString)
    field_id.SetWidth(14)
    newLayer.CreateField(field_id) 
    field_city_name = ogr.FieldDefn("city_name", ogr.OFTString)
    field_city_name.SetWidth(45)
    newLayer.CreateField(field_city_name)
    field_city_code = ogr.FieldDefn("city_code", ogr.OFTString)
    field_city_code.SetWidth(3)
    newLayer.CreateField(field_city_code)

    # Dereference the feature
    newFeature = None

# Dereference source
newDataSource = None

print(f"\nCreating .shp :\n")

# Step 3 : insert features in respectives cities layers

# List of layers create by city codes
os.chdir(outputFolder)
indexFile = 0
for i, file in enumerate(glob.glob("*.shp")):
    indexFile += 1
    currentDataSource = file
    currentLayer = file[0:-4]
    currentCodeCom = file[10:13]    

    # Open each communal layer in writing mode
    fileCommunal = f"{currentLayer}"
    inputFileCommunal = outputFolder + fileCommunal + extension
    dataSourceCommunal = driver.Open(inputFileCommunal, 1) # 0 read-only, 1 write
    layerCommunal = dataSourceCommunal.GetLayer()

    # Selection by city code in source layer
    sqlCommunal = f'SELECT * FROM {myFile} WHERE city_code LIKE \'{currentCodeCom}\''
    layerQueryCommunal = dataSource.ExecuteSQL(sqlCommunal)

    # Get attributes
    for feature in layerQueryCommunal :
        id = feature.GetField("id")        
        city_name = feature.GetField("city_name")
        city_code = feature.GetField("city_code")
        geom = feature.GetGeometryRef()

        # Create a new feature
        newFeatureCommunal = ogr.Feature(layerCommunal.GetLayerDefn())

        # Insert values
        newFeatureCommunal.SetField("id", id)        
        newFeatureCommunal.SetField("city_name", city_name)
        newFeatureCommunal.SetField("city_code", city_code)
        newFeatureCommunal.SetGeometry(geom)        

        # Insertfeature
        layerCommunal.CreateFeature(newFeatureCommunal)

        # Dereference the feature
        newFeatureCommunal = None

    # Dereference source
    dataSourceCommunal = None

    print(f"{indexFile}/{numberCityCode} - {file}")

print(f"\n{indexFile} layer created\n\n---- End of process ----")

Maybe this method may not be the best, in a Pythonic and syntactic point of view. Using functions could be a better way. If someone has better to offer, I am interested, especially in terms of performance and execution time (which is a big deal with layers containing a lot of features).

| improve this answer | |

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.