Michael,

You need to attach the CRS at layer creation time with:

layer = ds.CreateLayer(layer_name, geom_type=ogr.wkbPolygon, srs=sr)

Setting it with featureDefn.GetGeomFieldDefn(0).SetSpatialRef() is illegal. Now that https://gdal.org/development/rfc/rfc97_feature_and_fielddefn_sealing.html has been implemented, you should have got an error... (assuming you run your script with recent GDAL master)

from osgeo import ogr
ogr.UseExceptions()
ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource('/tmp/foo.shp')
lyr = ds.CreateLayer('foo')
from osgeo import osr
sr = osr.SpatialReference()
sr.SetFromUserInput("OGC:CRS84")
featureDefn = lyr.GetLayerDefn()
featureDefn.GetGeomFieldDefn(0).SetSpatialRef(sr)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/even/gdal/gdal/build_cmake/swig/python/osgeo/ogr.py", line 6017, in SetSpatialRef
    return _ogr.GeomFieldDefn_SetSpatialRef(self, *args)
RuntimeError: OGRGeomFieldDefn::SetSpatialRef() not allowed on a sealed object

I can't think of a reason why your script wouldn't work with a FlatGeoBuf output (did you get an error message?), which would be much better, as you could also attach the XRES, YRES, BANDCOUNT, etc. metadata items that would help the GTI driver avoid probing one of the tiles.

Even

Le 30/01/2024 à 12:31, Michael Sumner a écrit :
awesome, thanks Even I'm having fun with this one.

For anyone interested I created Python to parse the OpenTopography COP90 VRT (I have to wget it locally as I don't know how to hit the URL for the xml yet).

https://github.com/mdsumner/cog-example/blob/1ca74f3baffe69830180031ddabcbbc569816150/gti/cop90.py

(I'm using the DstRect to get efficient extent polygon of each tif without opening them).

This writes a cop90.shp which can be GTI'd via

gdalinfo GTI:/vsicurl/https://github.com/mdsumner/cog-example/raw/main/gti/cop90.shp

(seems I didn't get the SRS attached correctly).

I also have an R script in there for COP30 that I did previously to FlatGeoBuf, I couldn't get FlatGeoBuf to work from Python and don't know why.

If there's any pythonic-idiomatic advice on my code I'm very interested, slowly I'll try to use this to do the same in C++.  If there's a faster way to get these extents with GDAL I'm also keen to hear, thanks!

Cheers, Mike

On Fri, Jan 26, 2024 at 10:19 PM Even Rouault via gdal-dev <gdal-dev@lists.osgeo.org> wrote:

    Hi,

    the driver has recently landed into GDAL master, renamed as GTI:
    https://gdal.org/drivers/raster/gti.html

    Can be tested using
    https://gdal.org/download.html#gdal-master-conda-builds , the
    refreshed
    ghcr.io/osgeo/gdal <http://ghcr.io/osgeo/gdal> Docker image, or
    other nightly builds of GDAL master

    Even

    Le 20/12/2023 à 19:33, Even Rouault via gdal-dev a écrit :
    > Hi,
    >
    > For those not actively following github tickets & PR, I just
    want to
    > point to a new pending major functionality to improve management of
    > virtual mosaics with a very large number of tiles/sources (>
    tens of
    > thousands of tiles), by referencing them as features of a vector
    layer
    > (typically created by gdaltindex), instead of a XML file as in
    > traditional VRT, augmented with additional metadata.
    >
    > More details in https://github.com/OSGeo/gdal/pull/8983 (and in
    > initial ticket in https://github.com/OSGeo/gdal/issues/8861)
    >
    > Even
    >
-- http://www.spatialys.com
    My software is free, but my time generally not.

    _______________________________________________
    gdal-dev mailing list
    gdal-dev@lists.osgeo.org
    https://lists.osgeo.org/mailman/listinfo/gdal-dev



--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com

--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to