[Gdal-dev] IntegerList datatype

[Frank Warmerdam]

nida Khan wrote:
> Hi Everyone,
> 
> I'm in the process of coverting .NTF files to postgres table, ogr2ogr works
> great!
> 
> For one of the file, og2ogr creates a file with some messages like,
> 
> Cant create filed "Geom_id_link" with type IntegerList on PostgreSql layers,
> Creating as Varchar.
> 
> when I run sql on that table, I get this data.
> (3:150,340,2)
> (4,33,545,434,34)
> (10,232,354,23,343,...)
> 
> has anyone ever encountered this sort of problem?

Nida,

It is a fundamental flaw in the NTF readers in OGR that they don't assemble
the geometries for polygons from the component line strings.  I tried the
following python script using OGR on BL2000 data and it sort of worked, but
any of the long geom_id_of_link lists are truncated and ended with ... so the
data you have in the database now can never be properly assembled into
polygons.

The script can be run against a raw NTF file but in it's currently unfinished
state it does not actually write the assembled features anywhere, so it isn't
in a useful state yet.  Think of this as an experiment on my part in using
Scripting for processing features.

I will look at completing the script tonight, but no promises. Alternatively
I may just fix the NTF reader to do the same thing internally.

Best regards,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

import osr
import ogr
import string

#ds = ogr.Open( '/u/data/ntf/bl2000/HALTON.NTF' )
ds = ogr.Open( 'PG:dbname=test', update = 1 )

layer_count = ds.GetLayerCount()

#############################################################################-
# Establish access to the line and polygon layers.  Eventually we shouldn't
# hardcode this.

line_layer = ds.GetLayer(0)
poly_layer = ds.GetLayer(1)

#############################################################################
# Read all features in the line layer, holding just the geometry in a hash
# for fast lookup by GEOM_ID.

lines_hash = {}

feat = line_layer.GetNextFeature()
geom_id_field = feat.GetFieldIndex( 'GEOM_ID' )
tile_ref_field = feat.GetFieldIndex( 'TILE_REF' )
while feat is not None:
     geom_id = feat.GetField( geom_id_field )
     tile_ref = feat.GetField( tile_ref_field )

     if not lines_hash.has_key( tile_ref ):
         lines_hash[tile_ref] = {}

     sub_hash = lines_hash[tile_ref]
     sub_hash[geom_id] = feat.GetGeometryRef().Clone()

     feat.Destroy()

     feat = line_layer.GetNextFeature()

print 'Got %d lines.' % len(lines_hash)


#############################################################################
# Read all polygon features.

feat = poly_layer.GetNextFeature()
link_field = feat.GetFieldIndex( 'GEOM_ID_OF_LINK' )
tile_ref_field = feat.GetFieldIndex( 'TILE_REF' )

while feat is not None:
     tile_ref = feat.GetField( tile_ref_field )
     link_list = feat.GetField( link_field )

     # If the list is in string form we need to convert it.
     if type(link_list).__name__ == 'str':
         colon = string.find(link_list,':')
         items = string.split( link_list[colon+1:-1], ',' )
         link_list = []
         for item in items:
             try:
                 link_list.append(int(item))
             except:
                 print 'item failed to translate: ', item

     link_coll = ogr.Geometry( type = ogr.wkbGeometryCollection )
     for geom_id in link_list:
         geom = lines_hash[tile_ref][geom_id]
         link_coll.AddGeometry( geom )

     try:
         poly = ogr.BuildPolygonFromEdges( link_coll )
         print poly.ExportToWkt()
         feat.SetGeometryDirectly( poly )
     except:
         print 'BuildPolygonFromEdges failed.'

#    poly_layer.SetFeature( feat )
     feat.Destroy()

     feat = poly_layer.GetNextFeature()