[gdal-dev] GML and Spatial Reference

[Even Rouault]

Le vendredi 02 mai 2014 20:59:58, Ochs, Elke ERDC-RDE-CRREL-NH a écrit :
> Hello,
> 
> I'm not able to save a spatial reference to a GML file and would like help
> understanding why.  I'm creating a new GML from a copied feature.  I
> assign the spatial reference to both the output layer and geometry. 
> During debug, the layer and geometry report  that the SR has been assigned
> correctly (GetSratialRef/GetSpatialReference).  When I open the result GML
> file in another script, the spatial reference is empty.

Well, this is tricky :
- SRS in GML are only supported as EPSG codes. So you cannot pass an arbitrary 
proj.4 string
- Even if you do so, you will find that reading back the file will not lead to a 
SRS to be retrieved, since the driver is conservative by default (a GML file 
could have a mix of various SRS and we cannot trust the first one found in the 
file).

> 
> It's a non-standard SRS (albers, units feet) - not sure if this is an
> issue.
> 
> I'm using GDAL 1.10 compiled with Expat.  I'm fairly new to all of this. 
> Any help is appreciated.  I'm including most of the script, b/c I'm not
> sure where my error might be.  Some of it is taken from the  GDAL/OGR
> Cookbook.
> 
>     drv = ogr.GetDriverByName("GML")
>     ds = drv.Open(gmlin,0)
> 
>     srs = osr.SpatialReference()
>     proj4str = "+proj=aea +lat_1=29.5n +lat_2=45.5n +lat_0=23.0n
> +lon_0=96.0w +x_0=0.0 +y_0=0.0 +units=us-ft +datum=NAD83"
> srs.ImportFromProj4(proj4str)
> 
>     lyr = ds.GetLayer()
>     lyr.SetAttributeFilter("cat = " + cat)
>     outds = drv.CreateDataSource(gmlout,options = ["FORMAT=GML3"])
>     outlyr = outds.CreateLayer(outlyrname,srs,lyr.GetGeomType())
>     inLayerDefn = lyr.GetLayerDefn()
>     for i in range(0, inLayerDefn.GetFieldCount()):
>         fieldDefn = inLayerDefn.GetFieldDefn(i)
>         outlyr.CreateField(fieldDefn)
>     outLayerDefn = outlyr.GetLayerDefn()
> 
>     for feat in lyr:
>         geom = feat.GetGeometryRef()
>         outfeat = ogr.Feature(outLayerDefn)
>         for i in range(0, outLayerDefn.GetFieldCount()):
>             outfeat.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(),
> feat.GetField(i)) geom.AssignSpatialReference(srs)
>         outfeat.SetGeometry(geom)
>         outlyr.CreateFeature(outfeat)
> 
> Thank you,
> Elke

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html