Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I've got a similar question to this one How to programatically populate a shapefile field with polygon areas in OGR?, but the solution doesn't seem to be working for me. I am trying to import a shapefile from PostGIS into python using OGR, add a field, then populate the field with data from a csv based on an ID and the date. When I print the fields for the feature during the loop, the field gets the correct data but it isn't writing to the layer o3_proj (last two print statements). Am I using the SetFeature function incorrectly? I'm new to PostGIS and ogr, so please speak slowly!

import os, ogr, gdal, csv
driver = ogr.GetDriverByName('PostgreSQL')
datasource = driver.Open("PG: host='localhost' port='54321' dbname='chicago' user='postgres'", update = 1)

o3_proj = datasource.GetLayerByName('o3_monitors(geom3528)')

date = []
maxozone = []
aqi = []
aqs_site_id = []

f = csv.reader(open('/ozone/ozone_data.csv','rU'))
x = list(f)

for iii in range(1,len(x)):
    date.append(x[iii][0])
    aqs_site_id.append(x[iii][1])
    maxozone.append(float(x[iii][2]))
    aqi.append(int(x[iii][3]))

daily_field = ogr.FieldDefn('temp_o3',ogr.OFTReal)
o3_proj.CreateField(daily_field)

for feat in o3_proj:
    for jjj in range(0,len(date)):
        if aqs_site_id[jjj] == feat.GetField(0) and date[jjj] == '03/06/2005':
            o3_proj.SetFeature(feat)
            feat.SetField2(6,maxozone[jjj])
            o3_proj.SetFeature(feat)            
            print str(feat.GetField(0))+','+str(feat.GetField(1))+','+str(feat.GetField(2))+','+str(feat.GetField(3))+','+str(feat.GetField(4))+','+str(feat.GetField(5))+','+str(feat.GetField(6))
            print o3_proj.GetFeature(feat.GetFID()).GetField(6)

Result of the print statements:

17-031-0064,41.79079,-87.60165,2001,2012,None,0.036
None
17-031-0076,41.7514,-87.71349,2004,2012,None,0.038
None
17-031-4007,42.06028,-87.86323,2002,2012,None,0.035
None
17-031-4201,42.14,-87.79922,2001,2012,None,0.031
None
share|improve this question
    
Maybe you can simplify your code and variables names a little bit. Just enough to reproduce to problem once. But I think you are missing the CreateFeature statement. After you set a new field, you have to create the feature. Have a look at this example: pcjericks.github.io/py-gdalogr-cookbook/… –  ustroetz Apr 6 at 9:43
add comment

2 Answers

I suggest adding to the end of your script:-

o3_proj.Destroy()
o3_proj = None

As from other OGR tutorials I've read when working on my scripts, these calls at the end are suggested as necessary to ensure not just memory release, but also writing out data.

Let us know how it goes - as your use of SetFeature() does appear correct as far as I can tell. And I can see you've opened the file in update mode.

The other possibility I can think of, if my suggestion doesn't work - perhaps updating isn't supported for PostGIS files by the GDAL OGR driver?? You'd need to check.

share|improve this answer
add comment

I know this isn't a straight answer to your question but I advice you to use pyshp to create and manipulate shape files. It is written in pure python and very easy to use. To edit a shapefile simply do:

e = shapefile.Editor(shapefile="shapefiles/test/name.shp")

Then add a field

e.field("add_field","C","30")

Well, then kick off your loop and populate it by

e.record("value")

When you are done, call the save method w.save(). That's it.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.