Processing feature string fields to raster

Question

I am using the GDAL python bindings. I am trying to rasterize a layer from an arcinfo coverage file, but the values I need to rasterize I have to calculate from string fields. OGR can not write to arcinfo coverage files and, in any event, I do not want to save the newly calculated fields except to the raster. Would the following work?

• Load coverage layer
• Create OGR dataset of type memory
• Copy layer to memory dataset
• Iterate over features in memory dataset, calculating new field and adding to feature
• Rasterize memory dataset

Or is there a better (or more simple) method I am overlooking?

Answer 1

After a day of struggling with GDAL, I can confirm the methodology works and runs quickly. I was working with ice chart E00 files from CIS which I had converted to arcinfo coverage using avcimport. I have posted my prototype code below to help others. The E00 was converted to an arcinfo coverage called "test". The script must be run from the same directory i.e. the parent directory to "test". I am new to python and GIS and the code is a bit ugly and missing error checking, but at least shows how it can be done. I wanted ice concentration and an estimate of ice thickness from the ice stages. Land fast ice and land were assigned values of -1 and -2. The coordinate transformation was a simple rotation of the Lambert Conformal Projection to fit the data more nicely into the raster i.e. minimal NoData pixels.

Code guaranteed only to work on ice chart for Hudson Bay Regional 1997年01月01日. Should also work for later dates, but untested and some charts have glitches.

from osgeo import gdal, ogr, osr
c_conc = { '': 0.0, '0.': 0.0, '1': 0.1, '2': 0.2, '3': 0.3, '4': 0.4, '5': 0.5,
 '6': 0.6, '7': 0.7, '8': 0.8, '9': 0.9, '9+': 1.0, '10': 1.0}
#s_thick = { '': 0.0, '1': 0.05, '2': 0.05, '4': 0.125, '5': 0.225, '6': 1.0,
# '7': 0.5, '1.': 0.95, '4.': 1.6, '7.': 2.0, '8.':2.0, '9.': 2.0}
s_thick = { '': 0.0, '1': 0.05, '4': 0.125, '5': 0.225, '7': 0.5,
 '1.': 0.95, '4.': 1.6, '7.': 2.0, '8.':2.0, '9.': 2.0}
# legend definitions 1 = use values , 2 = flag as land value, 3 = flag as Fast Ice
feature_legend = {'remote egg': 1, 'no data': 1, 'bergy water': 1, 'land': 2, 'fast ice': 3}
land_value = -2
fast_value = -1
pixel_size = 10000
avcbin_path = 'test'
raster_path = 'temp.tif'
drv = ogr.GetDriverByName('AVCBin')
ds = drv.Open(avcbin_path, 0)
layer = ds.GetLayerByName('PAL')
s_srs = layer.GetSpatialRef()
t_srs = s_srs.Clone()
err = t_srs.SetProjParm('central_meridian', -59.1)
srs_trans = osr.CoordinateTransformation(s_srs, t_srs)
feature = layer.GetNextFeature()
datestring = feature.GetFieldAsString('DATE_CARTE').strip()
if datestring == '0':
 print 'WARNING: bad datestring'
print avcbin_path, datestring
m_drv = ogr.GetDriverByName('Memory')
m_ds = m_drv.CreateDataSource('wrk')
m_ds.CopyLayer(layer, 'PAL')
ds = None
drv = None
layer = m_ds.GetLayerByName('PAL')
layer_defn = layer.GetLayerDefn()
conc_field = ogr.FieldDefn('CONC', ogr.OFTInteger)
thick_field = ogr.FieldDefn('THICK', ogr.OFTInteger)
layer.CreateField(conc_field)
layer.CreateField(thick_field)
#for feature in layer:
# print feature.GetFID(), feature.GetFieldAsString('DATE_CARTE'), \
# feature.GetFieldAsString('SOURCE'), feature.GetFieldAsString('A_LEGEND')
#m_ds = None
#m_drv = None
for feature in layer:
 geom = feature.GetGeometryRef()
 geom.Transform(srs_trans)
 datestring_tmp = feature.GetFieldAsString("DATE_CARTE").strip()
 egg_attr = feature.GetFieldAsString("EGG_ATTR")
 feature_key = feature.GetFieldAsString('A_LEGEND').strip().lower()
 if datestring_tmp != datestring:
 print feature.GetFID(), feature_key, egg_attr, datestring_tmp, datestring
 assert datestring_tmp == datestring
 try:
 feature_type = feature_legend[feature_key]
 ct = c_conc[feature.GetFieldAsString('E_CT').strip()]
 ca = c_conc[feature.GetFieldAsString('E_CA').strip()]
 cb = c_conc[feature.GetFieldAsString('E_CB').strip()]
 cc = c_conc[feature.GetFieldAsString('E_CC').strip()]
 cd = c_conc[feature.GetFieldAsString('E_CD').strip()]
 sa = s_thick[feature.GetFieldAsString('E_SA').strip()]
 sb = s_thick[feature.GetFieldAsString('E_SB').strip()]
 sc = s_thick[feature.GetFieldAsString('E_SC').strip()]
 sd = s_thick[feature.GetFieldAsString('E_SD').strip()]
 except:
 print feature_key, egg_attr
 m_ds = None
 m_drv = None
 raise
# finally:
# m_ds = None
# m_drv = None
 assert cd == 0.
 if feature_type == 1:
 if ca == 0.:
 ca = ct
 if ct == 0. and ca != 0.:
 ct = ca
 if abs(ct - (ca + cb + cc + cd)) > 0.01 and sd != 0.:
 cd = ct - (ca + cb + cc)
 concentration = int(ct * 10)
 if ct > 0.:
 thickness = int((ca*sa + cb*sb + cc*sc + cd*sd) / ct * 100)
 else:
 thickness = 0
 elif feature_type == 2:
 concentration = land_value
 thickness = land_value
 elif feature_type == 3:
 concentration = fast_value
 thickness = fast_value
 print feature.GetFID(), feature_key, feature_type, egg_attr, concentration, thickness
 feature.SetField('CONC', concentration)
 feature.SetField('THICK', thickness)
 layer.SetFeature(feature)
print 'Done' 
#layer.ResetReading()
#for feature in layer:
# print feature.GetFID(), feature.GetFieldAsString('A_LEGEND').strip().lower(), \
# feature.GetFieldAsInteger('CONC'), feature.GetFieldAsInteger('THICK')
# Rasterize
x_min, x_max, y_min, y_max = layer.GetExtent()
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
t_ds = gdal.GetDriverByName('GTiff').Create(raster_path, x_res, y_res, 2, gdal.GDT_Int16)
t_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
conc_band = t_ds.GetRasterBand(1)
conc_band.SetNoDataValue(-32768)
thick_band = t_ds.GetRasterBand(2)
thick_band.SetNoDataValue(-32768)
gdal.RasterizeLayer(t_ds, [1], layer, options=['ATTRIBUTE=CONC'])
gdal.RasterizeLayer(t_ds, [2], layer, options=['ATTRIBUTE=THICK'])
t_ds = None
m_ds = None
m_drv = None