Re: [gdal-dev] can GDAL render two images as if over-plotted?

[Michael Sumner via gdal-dev]

Just an update, I have found out how to do this numerically.Will
consider a PR for gdaldem if I can make it work at some point.

Thanks, Mike



On Tue, Jul 9, 2024 at 7:31 AM Michael Sumner  wrote:

> Is it possible to "layer" one image over another, the first at full
> transparency and the second with partial transparency?
>
> I'm thinking of combining a hillshade and a (partially) transparent
> color-relief image, a visual enhancement technique.   (I know how to do
> this downstream in R and Python or similar, I'd like to stay in C++'s warm
> embrace. )
>
> I can go and find out what "over plot" means when it comes to colour
> imagery, just wondering if there's anything in-built already.   This might
> be a nice additional mode to gdaldem (and is what I thought "color-relief"
> meant,  at first).
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Research Software Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>


-- 
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] can GDAL render two images as if over-plotted?

[Michael Sumner via gdal-dev]

Is it possible to "layer" one image over another, the first at full
transparency and the second with partial transparency?

I'm thinking of combining a hillshade and a (partially) transparent
color-relief image, a visual enhancement technique.   (I know how to do
this downstream in R and Python or similar, I'd like to stay in C++'s warm
embrace. )

I can go and find out what "over plot" means when it comes to colour
imagery, just wondering if there's anything in-built already.   This might
be a nice additional mode to gdaldem (and is what I thought "color-relief"
meant,  at first).

Cheers, Mike


-- 
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalwarp across the antimeridian with a source extent >180E

[Michael Sumner via gdal-dev]

Thanks!   Educational as ever, this closes out a few hanging threads for
me, much appreciated.

Cheers, Mike



On Mon, Jun 3, 2024 at 3:48 AM Even Rouault 
wrote:

> Michael,
>
> https://github.com/OSGeo/gdal/pull/10108 will fix it.
>
> You can also workaround the issue by overriding the extent of the source
> dataset to be exactly -179.995,89.995,180.005,-89.995 with
>
>
>  "vrt://NETCDF:2024010209-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:
> analysed_sst?a_srs=EPSG:4326_ullr=-179.995,89.995,180.005,-89.995"
>
> The issue was twofolds:
>
> - the dataset uses lon/lat arrays with single precision floating-point
> numbers. Hence -179.99 gets actually read a a   read as -179.99000549316406
> as a double-precision number, and thus the maxx - minxx difference was
> slightly over 360 degrees. I've added a heuristics to try to round
> -179.99-single-precision as -179.99-double-precision
>
> - which defeated a heuristics in the warper to identify the center
> longitude of a dataset registered in a geographic CRS. This center
> longitude is just (min_lon + max_lon) / 2 if the dataset doesn't cover more
> than 360 degrees of longitude. This value, when computed, is passed as a
> hint OGRCoordinateTransformation so that it can post-correct longitudes to
> apply a +/- 360 degree offset, to be in the range of the source dataset.
> I've relaxed the sanity check to allow slighly more than 360 degrees
>
> Even
>
>
> Le 01/06/2024 à 08:19, Michael Sumner via gdal-dev a écrit :
>
> This data source has an odd georeferencing, it's a 36000x17999 raster in
> the extent -179.995 -89.995 180.005 89.995.
>
> vrt://NETCDF:/vsicurl/
> https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/2024010109-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326
>
> (Why is it 17999: well there's one missing row - you'd expect 18000 at
> 0.01 resolution and that manifests as half a cell short at the top and
> bottom).
>
> But also the x registration is half a cell to the east from what you would
> expect.
>
> It's correctly registered in that frame though, if we shift it west by the
> half cell I've verified visually that it "looks wrong" compared to
> coastlines etc.
>
> The point of this report is, gdalwarp misses the western window of data
> for a request across the antimeridian (there's no data in the output east
> of 180).
>
> gdalwarp $dsn -te 173.00 -19.85 185.00 -15.00 out1.tif -co COMPRESS=LZW
>
> If we "fix" the extent to be more longitudinally pleasing, we get the data
> east of the antimeridian, the map is properly filled either side of the
> antimeridian.
>
>
> vrt://NETCDF:/vsicurl/
> https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/2024010109-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326_ullr=-180,89.995,180,-89.995
>
> gdalwarp $dsnfix -te 173.00 -19.85 185.00 -15.00 out2.tif -co COMPRESS=LZW
>
>
> Which is good, but as mentioned above the workaround means we are now half
> a cell out. It works fine for projected windows over the antimeridian in
> either situation.
>
> The compressed file sizes reflect the missing data in the second one,
> these are 625K and 1.1Mb.
>
> I'd like the warper to be able to correctly fill the data from the western
> window in the original case, we have determined that we can't feasibly
> "fix" the extent.
>
> Or maybe something else I'm missing? I don't *think* SAMPLE_GRID or
> SOURCE_EXTRA helps here. I know that I could craft a meridian-crossing
> combination in VRT but I'd rather like this workflow to work as-is.
> (Chasing down the weird grid referencing is another project, but it's a
> really massive dataset so I'm pretty unconfident of that occurring).
>
> Cheers, Mike
>
> --
> Michael Sumner
> Research Software Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalwarp across the antimeridian with a source extent >180E

[Michael Sumner via gdal-dev]

I forgot to mention that it needs earthdata credentials set up, basically
your "Authorization: Bearer " in GDAL_HTTP_HEADERS or similar
config.

https://urs.earthdata.nasa.gov/documentation/for_users/user_token

You can't download or stream these files without that set or logging in
(the file used is indexed here:
https://cmr.earthdata.nasa.gov/virtual-directory/collections/C1996881146-POCLOUD/temporal/2024/01/01
)

I also crafted a gist to keep the formatting easier to copy paste:

https://gist.github.com/mdsumner/4ab86d52dd981621ad508297ba1a409c

Cheers, Mike



On Sat, Jun 1, 2024 at 4:19 PM Michael Sumner  wrote:

> This data source has an odd georeferencing, it's a 36000x17999 raster in
> the extent -179.995 -89.995 180.005 89.995.
>
> vrt://NETCDF:/vsicurl/
> https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/2024010109-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326
>
> (Why is it 17999: well there's one missing row - you'd expect 18000 at
> 0.01 resolution and that manifests as half a cell short at the top and
> bottom).
>
> But also the x registration is half a cell to the east from what you would
> expect.
>
> It's correctly registered in that frame though, if we shift it west by the
> half cell I've verified visually that it "looks wrong" compared to
> coastlines etc.
>
> The point of this report is, gdalwarp misses the western window of data
> for a request across the antimeridian (there's no data in the output east
> of 180).
>
> gdalwarp $dsn -te 173.00 -19.85 185.00 -15.00 out1.tif -co COMPRESS=LZW
>
> If we "fix" the extent to be more longitudinally pleasing, we get the data
> east of the antimeridian, the map is properly filled either side of the
> antimeridian.
>
>
> vrt://NETCDF:/vsicurl/
> https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/2024010109-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326_ullr=-180,89.995,180,-89.995
>
> gdalwarp $dsnfix -te 173.00 -19.85 185.00 -15.00 out2.tif -co COMPRESS=LZW
>
>
> Which is good, but as mentioned above the workaround means we are now half
> a cell out. It works fine for projected windows over the antimeridian in
> either situation.
>
> The compressed file sizes reflect the missing data in the second one,
> these are 625K and 1.1Mb.
>
> I'd like the warper to be able to correctly fill the data from the western
> window in the original case, we have determined that we can't feasibly
> "fix" the extent.
>
> Or maybe something else I'm missing? I don't *think* SAMPLE_GRID or
> SOURCE_EXTRA helps here. I know that I could craft a meridian-crossing
> combination in VRT but I'd rather like this workflow to work as-is.
> (Chasing down the weird grid referencing is another project, but it's a
> really massive dataset so I'm pretty unconfident of that occurring).
>
> Cheers, Mike
>
> --
> Michael Sumner
> Research Software Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>


-- 
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] gdalwarp across the antimeridian with a source extent >180E

[Michael Sumner via gdal-dev]

This data source has an odd georeferencing, it's a 36000x17999 raster in
the extent -179.995 -89.995 180.005 89.995.

vrt://NETCDF:/vsicurl/
https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/2024010109-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326

(Why is it 17999: well there's one missing row - you'd expect 18000 at 0.01
resolution and that manifests as half a cell short at the top and bottom).

But also the x registration is half a cell to the east from what you would
expect.

It's correctly registered in that frame though, if we shift it west by the
half cell I've verified visually that it "looks wrong" compared to
coastlines etc.

The point of this report is, gdalwarp misses the western window of data for
a request across the antimeridian (there's no data in the output east of
180).

gdalwarp $dsn -te 173.00 -19.85 185.00 -15.00 out1.tif -co COMPRESS=LZW

If we "fix" the extent to be more longitudinally pleasing, we get the data
east of the antimeridian, the map is properly filled either side of the
antimeridian.


vrt://NETCDF:/vsicurl/
https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/2024010109-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326_ullr=-180,89.995,180,-89.995

gdalwarp $dsnfix -te 173.00 -19.85 185.00 -15.00 out2.tif -co COMPRESS=LZW


Which is good, but as mentioned above the workaround means we are now half
a cell out. It works fine for projected windows over the antimeridian in
either situation.

The compressed file sizes reflect the missing data in the second one, these
are 625K and 1.1Mb.

I'd like the warper to be able to correctly fill the data from the western
window in the original case, we have determined that we can't feasibly
"fix" the extent.

Or maybe something else I'm missing? I don't *think* SAMPLE_GRID or
SOURCE_EXTRA helps here. I know that I could craft a meridian-crossing
combination in VRT but I'd rather like this workflow to work as-is.
(Chasing down the weird grid referencing is another project, but it's a
really massive dataset so I'm pretty unconfident of that occurring).

Cheers, Mike

-- 
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] clear band level metadata

[Michael Sumner via gdal-dev]

Hello, I have this process to convert a URL to a MEM dataset.

It's not the target data but it's representative, a netcdf with band-level
and dataset-level metadata. We can clear the dataset level with
COPY_SRC_MDD, but I can't see how to do that for the band level. We're
hoping to keep this  without generating files, the final step uses VSI to
read into a byte array.

Can I clear the band metadata without an intermediate dataset? Opening in
update mode breaks the COG layout (or is that ok to do this with
IGNORE_COG_LAYOUT_BREAK=YES ?).


from osgeo import gdal
gdal.UseExceptions()

dsn = "vrt:///vsicurl/
https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202405/oisst-avhrr-v02r01.20240511.nc?sd_name=sst
"

opts = gdal.TranslateOptions(format = "COG",
 creationOptions = [ "COMPRESS=ZSTD", "PREDICTOR=STANDARD",
"RESAMPLING=AVERAGE",
 "SPARSE_OK=YES", "COPY_SRC_MDD=NO"])

ds = gdal.Open(dsn)
tmp = "/vsimem/tmp1.tif"
## here we lose the dataset-level metadata as desired with COPY_SRC_MDD
gdal.Translate(tmp, ds, options = opts)

## but the band metadata is still there
lds = gdal.Open(tmp)
#{'add_offset': '0', 'long_name': 'Daily sea surface temperature',
'NETCDF_DIM_time': '16932', 'NETCDF_DIM_zlev': '0', 'NETCDF_VARNAME':
'sst', 'scale_factor': '0.009998', 'units': 'Celsius', 'valid_max':
'4500', 'valid_min': '-300', '_FillValue': '-999'}
lds.Close()
ds.Close()



I guess I need to open and translate to another COG dataset after breaking
the COG layout?



## now open in rw, to zap band level md (but breaks COG layout)
#lds = gdal.OpenEx(tmp, gdal.GA_Update, open_options =
 ["IGNORE_COG_LAYOUT_BREAK=YES"])
#lds.GetRasterBand(1).SetMetadata({})
#lds.Close()
#Warning 1: tmp1.tif: The IFD has been rewritten at the end of the file,
which breaks COG layout.

ds = gdal.Open(tmp)
#Warning 1: tmp1.tif: This file used to have optimizations in its layout,
but those have been, at least partly, invalidated by later changes

Cheers, Mike




--
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] another problem HDF5

[Michael Sumner via gdal-dev]

Thanks Even!I've reported it to the providers.

(btw don't gazump me on that book)

Cheers, Mike


On Tue, Apr 30, 2024 at 11:04 PM Even Rouault 
wrote:

> Hi,
>
> oh well, so this i actually a HDF5 file using netCDF conventions, but
> likely created with the HDF5 API itself. Per se, this isn't the issue, as
> you've figured it, we can convince the netCDF driver to open it. The lack
> of CRS comes from the fact that the land_mask_map variable has this
> attribute:
>
> land_mask_map:grid_mapping = "crs: grid_x crs: grid_y" ;
>
> Which tries to point to the "crs" variable where the CRS is defined. In
> 99% of the products, the value of grid_mapping would be just "crs", and
> life would be good.
>
> Here they use the "second format" for grid_mapping, described at
> https://cfconventions.org/Data/cf-conventions/cf-conventions-1.11/cf-conventions.html#grid-mappings-and-projections
>
> "In the second format, it is a blank-separated list of words
> ":  [ …​]
> [: …​]", which identifies one or
> more grid mapping variables, and with each grid mapping associates one or
> more coordinatesVariables, i.e. coordinate variables or auxiliary
> coordinate variables."
>
> That documentation isn't 100% crystal clear but the later example 5.10
> "British National Grid" gives an example:
>
> temp:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
>
> So you can point to several crs variables, but for each you need to define
> *all* the relevant axis. The way it is done in the product you point
> doesn't match that philosophy (if they wanted to use the second format that
> should rather set "crs: grid_x grid_y", but using the second format here is
> just asking for troubles), and GDAL doesn't manage to parse it.
>
> You could tune netCDFDataset::SetProjectionFromVar() to recognize this
> weird value for grid_mapping in netCDFDataset::SetProjectionFromVar(),
> around line 3287, but I don't feel myself personally motivated to do that,
> if that product is just an outlier, so I leave it to you as an exercise if
> you wish :-). I guess just testing against the hardcoded "crs: grid_x crs:
> grid_y" value would be good enough to handle that outlier.
>
> (Someone, super bored, should author a book about with horror stories with
> netCDF and HDF georeferencing. There's a lot of material. Although likely
> not to be a best seller)
>
> Even
> Le 30/04/2024 à 09:45, Michael Sumner via gdal-dev a écrit :
>
> This time, with:
>
>
> https://n5eil01u.ecs.nsidc.org/ATLAS/ATL20.004/2018.10.14/ATL20-01_20181001010615_00370101_004_01.h5
>
> NETCDF gets the geotransform (from x_grid,y_grid which report in the
> metadata) but no crs:
>
>   gdalinfo NETCDF:ATL20-02_20181001010615_00370101_004_01.h5 -sd 1 -nomd
> Driver: netCDF/Network Common Data Format
> Files: ATL20-02_20181001010615_00370101_004_01.h5
> Size is 316, 332
> Origin = (-395.000,435.000)
> Pixel Size = (25000.000,-25000.000)
> Corner Coordinates:
> Upper Left  (-395.000, 435.000)
> Lower Left  (-395.000,-395.000)
> Upper Right ( 395.000, 435.000)
> Lower Right ( 395.000,-395.000)
> Center  (   0.000,  20.000)
> Band 1 Block=500x500 Type=Float64, ColorInterp=Undefined
>   Unit Type: degrees_north
>
>
> HDF5 doesn't get geotransform or crs:
>
> gdalinfo ATL20-02_20181001010615_00370101_004_01.h5 -sd 1 -nomd
> Driver: HDF5Image/HDF5 Dataset
> Files: ATL20-02_20181001010615_00370101_004_01.h5
> Size is 316, 332
> Corner Coordinates:
> Upper Left  (0.0,0.0)
> Lower Left  (0.0,  332.0)
> Upper Right (  316.0,0.0)
> Lower Right (  316.0,  332.0)
> Center  (  158.0,  166.0)
> Band 1 Block=500x500 Type=Float32, ColorInterp=Undefined
>   NoData Value=3.4028235e+38
>
> This is unproblematically NSIDC 25km south, there's no mix of poles in
> this one. I honestly don't know what to expect from these files generally
> though I get my hopes up sometimes haha.
>
> I'll have to keep notes and explore more deeply with native HDF5 tools.
>
> Cheers, Mike
>
>
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdal2tiles for floating point images

[Michael Sumner via gdal-dev]

Ok so my naive edits are clearly not enough, they still get written as Byte
so it's deeper in the target spec and spread across a few places I'm not
ready to get across yet.

Happy to pursue in the longer term though.

Cheers, Mike


On Tue, Apr 30, 2024 at 7:00 AM Michael Sumner  wrote:

>
>
> On Fri, Apr 26, 2024 at 5:37 AM lefsky--- via gdal-dev <
> gdal-dev@lists.osgeo.org> wrote:
>
>> I'd like to have a version of gdal2tiles that handles image types other
>> than uint8. Is there a reason why it doesn't handle those types of images?
>> It appears I'd have to modify the output format from png to tiff.  I have
>> never modified gdal source before and I'm wondering if this modification
>> (which appears straightforward to a naive user) hasn't been done before.
>>
>
>
> Hello, it appears to be sufficient to add GTiff to tiledriver options,
> handle the file extension, and add an option to override the error/message
> that advises conversion to Byte:
>
>
> https://github.com/OSGeo/gdal/compare/master...dis-organization:gdal:gdal2tiles-nonbyte
>
> I wanted this myself as I'm working on a version of the tiling in an R
> package, and I'm not confident enough about it without being able to
> compare to gdal2tiles.py.
>
> It will need care to prepare it and merge into GDAL itself, with some
> thought about documentation, implications for html output (it's compleltely
> incompatible with the html produced I expect), and tests, but if you want
> to use it locally it seems ok and I'm happy to follow up off-list.
>
> Note that you can use VRT to reference such a tiled output for numeric
> data, I used it for a tiled global elevation source here (so if we fold
> this into GDAL that could be a nice option to include, VRT for reading as
> TMS):
>
> https://github.com/hypertidy/sds/blob/main/R/sources.R#L41
>
> Cheers, Mike
>
>
>
>
>>
>> Michael
>>
>>
>> --
>> Michael Lefsky (He/His)
>> Home Location: HVHF+GH
>> Cell: 970-980-9036
>> http://www.researcherid.com/rid/A-7224-2009
>>
>> *“for being prematurely, and worse, intuitively right — there’s a heavy
>> price. But for being wrong — no, not so long as you’re wrong in a pack."
>> Gary Brecher / Portis*
>>
>> *I acknowledge that I live and work on stolen land. This is the land of
>> the Cheyenne, Arapaho, Ute, and Ocheithi Sakowin people. To learn more
>> about these nations, please visit;
>> http://www.utemountainutetribe.com/
>> http://www.cheyennenation.com/
>> https://cheyenneandarapaho-nsn.gov/
>> https://native-land.ca/
>>
>> ___________
>> 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
>


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] another problem HDF5

[Michael Sumner via gdal-dev]

This time, with:

https://n5eil01u.ecs.nsidc.org/ATLAS/ATL20.004/2018.10.14/ATL20-01_20181001010615_00370101_004_01.h5

NETCDF gets the geotransform (from x_grid,y_grid which report in the
metadata) but no crs:

  gdalinfo NETCDF:ATL20-02_20181001010615_00370101_004_01.h5 -sd 1 -nomd
Driver: netCDF/Network Common Data Format
Files: ATL20-02_20181001010615_00370101_004_01.h5
Size is 316, 332
Origin = (-395.000,435.000)
Pixel Size = (25000.000,-25000.000)
Corner Coordinates:
Upper Left  (-395.000, 435.000)
Lower Left  (-395.000,-395.000)
Upper Right ( 395.000, 435.000)
Lower Right ( 395.000,-395.000)
Center  (   0.000,  20.000)
Band 1 Block=500x500 Type=Float64, ColorInterp=Undefined
  Unit Type: degrees_north


HDF5 doesn't get geotransform or crs:

gdalinfo ATL20-02_20181001010615_00370101_004_01.h5 -sd 1 -nomd
Driver: HDF5Image/HDF5 Dataset
Files: ATL20-02_20181001010615_00370101_004_01.h5
Size is 316, 332
Corner Coordinates:
Upper Left  (0.0,0.0)
Lower Left  (0.0,  332.0)
Upper Right (  316.0,0.0)
Lower Right (  316.0,  332.0)
Center  (  158.0,  166.0)
Band 1 Block=500x500 Type=Float32, ColorInterp=Undefined
  NoData Value=3.4028235e+38

This is unproblematically NSIDC 25km south, there's no mix of poles in this
one. I honestly don't know what to expect from these files generally though
I get my hopes up sometimes haha.

I'll have to keep notes and explore more deeply with native HDF5 tools.

Cheers, Mike




-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] HDF5 and geolocation arrays

[Michael Sumner via gdal-dev]

you need

gdalinfo --config GDAL_HTTP_HEADERS="Authorization: Bearer "
 "/vsicurl/https://n5eil01u.ecs.nsidc.org/…

Where "" is an Earthdata access token from your account.

https://urs.earthdata.nasa.gov/documentation/for_users/user_token

HTH   (there are other methods like setting GDAL_HTTP_HEADER_FILE to an
actual file with that "Authorization: Bearer ..." content.



On Tue, Apr 30, 2024 at 7:37 AM Joaquim Manuel Freire Luís 
wrote:

> > This HDF5 (requires earthdata credentials your "Authorization: Bearer
> " in GDAL_HTTP_HEADERS, or equiv) presents without geolocation
> arrays.
>
>
>
> gdalinfo "/vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5;
> -sd 26
>
>
>
> Excuse for this little derail, but how do we do that? I mean, the
> credentials. I tried with both:
>
>
>
> gdalinfo --config GDAL_HTTP_HEADERS=login:passw "/vsicurl/
> https://n5eil01u.ecs.nsidc.org/…
>
>
>
> and
>
>
>
> gdalinfo  "/vsicurl/https://login:pas...@n5eil01u.ecs.nsidc.org…
>
>
>
> but both errored with
>
>
>
> ERROR 3: Cannot read 4029 bytes
>
> gdalinfo failed - unable to open '/vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
> '.
>
>
>
> and
>
>
>
> ERROR 3: Cannot read 4029 bytes
>
> gdalinfo failed - unable to open '/vsicurl/
> https://login:pas...@n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
> '.
>
>
>
>
>
> Thanks
>
>
>
> Joaquim
>
>
>
> *From:* gdal-dev  *On Behalf Of *Michael
> Sumner via gdal-dev
> *Sent:* Monday, April 29, 2024 8:11 PM
> *To:* gdal-dev 
> *Subject:* [gdal-dev] HDF5 and geolocation arrays
>
>
>
> This HDF5 (requires earthdata credentials your "Authorization: Bearer
> " in GDAL_HTTP_HEADERS, or equiv) presents without geolocation
> arrays.
>
>
>
> gdalinfo "/vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5;
> -sd 26
> Driver: HDF5Image/HDF5 Dataset
> Files: /vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
> Size is 608, 896
> Metadata:
>   Conventions=CF-1.6
>   HDFEOS_INFORMATION_HDFEOSVersion=HDFEOS_5.1.15
>   history=This version of the Sea Ice processing code contains updates
> provided by the science team on September 16, 2019. For details on these
> updates, see the release notes provided in the DAP.
>   institution=NASA's AMSR Science Investigator-led Processing System (SIPS)
>   references=Please cite these data as: Markus, T., J. C. Comiso, and W.
> N. Meier. 2018. AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness
> Temperatures, Sea Ice Concentration, Motion & Snow Depth Polar Grids,
> Version 1. [Indicate subset used]. Boulder, Colorado USA. NASA National
> Snow and Ice Data Center Distributed Active Archive Center. doi:
> https://doi.org/10.5067/RA1MIJOYPK3P.
>   source=satellite observation
>   title=AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures, Sea
> Ice Concentration, Motion & Snow Depth Polar Grids
> Corner Coordinates:
> Upper Left  (0.0,0.0)
> Lower Left  (0.0,  896.0)
> Upper Right (  608.0,0.0)
> Lower Right (  608.0,  896.0)
> Center  (  304.0,  448.0)
> Band 1 Block=608x1 Type=Int32, ColorInterp=Undefined
>   Metadata:
> comment=data value meaning: 0 -- Open Water, 110 -- missing/not
> calculated, 120 -- Land
> coordinates=lon lat
> long_name=Sea ice concentration daily average
> units=percent
>
>
>
>
>
>
>
> gdalinfo --version
> GDAL 3.9.0dev-cb4d30f56d, released 2024/04/15
>
>
>
> The geolocation arrays are sds 33 and 32 respectively:
>
>
>
> HDF5:"/vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
> "://HDFEOS/GRIDS/NpPolarGrid12km/lon
>
>
>
> HDF5:"/vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
> "://HDFEOS/GRIDS/NpPolarGrid12km/lat
>
>
>
> And things work when lining those up in VRT with warp. Can the HDF5 driver
> be made to auto-detect these geolocation arrays?
>
>
>
> I see that the NETCDF driver actually does:
>
>
>
> gdalinfo "NetCDF:/vsicurl/
> https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5;
> -sd 26
>
>
>
> I'm asking as an email rather than pursuing the fix because, these data
> are actual

Re: [gdal-dev] gdal2tiles for floating point images

[Michael Sumner via gdal-dev]

On Fri, Apr 26, 2024 at 5:37 AM lefsky--- via gdal-dev <
gdal-dev@lists.osgeo.org> wrote:

> I'd like to have a version of gdal2tiles that handles image types other
> than uint8. Is there a reason why it doesn't handle those types of images?
> It appears I'd have to modify the output format from png to tiff.  I have
> never modified gdal source before and I'm wondering if this modification
> (which appears straightforward to a naive user) hasn't been done before.
>


Hello, it appears to be sufficient to add GTiff to tiledriver options,
handle the file extension, and add an option to override the error/message
that advises conversion to Byte:

https://github.com/OSGeo/gdal/compare/master...dis-organization:gdal:gdal2tiles-nonbyte

I wanted this myself as I'm working on a version of the tiling in an R
package, and I'm not confident enough about it without being able to
compare to gdal2tiles.py.

It will need care to prepare it and merge into GDAL itself, with some
thought about documentation, implications for html output (it's compleltely
incompatible with the html produced I expect), and tests, but if you want
to use it locally it seems ok and I'm happy to follow up off-list.

Note that you can use VRT to reference such a tiled output for numeric
data, I used it for a tiled global elevation source here (so if we fold
this into GDAL that could be a nice option to include, VRT for reading as
TMS):

https://github.com/hypertidy/sds/blob/main/R/sources.R#L41

Cheers, Mike




>
> Michael
>
>
> --
> Michael Lefsky (He/His)
> Home Location: HVHF+GH
> Cell: 970-980-9036
> http://www.researcherid.com/rid/A-7224-2009
>
> *“for being prematurely, and worse, intuitively right — there’s a heavy
> price. But for being wrong — no, not so long as you’re wrong in a pack."
> Gary Brecher / Portis*
>
> *I acknowledge that I live and work on stolen land. This is the land of
> the Cheyenne, Arapaho, Ute, and Ocheithi Sakowin people. To learn more
> about these nations, please visit;
> http://www.utemountainutetribe.com/
> http://www.cheyennenation.com/
> https://cheyenneandarapaho-nsn.gov/
> https://native-land.ca/
>
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] HDF5 and geolocation arrays

[Michael Sumner via gdal-dev]

Oh wow, that's fantastic - thanks so much. Lots for me to learn from there
too, with the  HDF5EOS GRID details.

Cheers, Mike

On Tue, Apr 30, 2024 at 6:08 AM Even Rouault 
wrote:

> Michael,
>
> actually support for those products was already there at 95% since they
> use the HDF5EOS GRID formalism which we support since 3.7. But that
> particular type of products has an oddity. They have an empty
> GROUP=Dimension within the GROUP=GridStructure in the HDF5EOS metadata,
> which our parser didn't like as we use the information in this group to map
> dimension names to dimension index and size. Fixed/worked around in
> https://github.com/OSGeo/gdal/pull/9807 . Hard to tell if those products
> are in fault given that the HDF5EOS specification at
> https://www.earthdata.nasa.gov/s3fs-public/imported/ESDS-RFC-008-v1.1.pdf
> , appendix A, is extremely minimum regarding the HDF5EOS metadata itself
> (for once, I complain for a spec to be too minimal !), and the example they
> give for grids doesn't make sense to me (the declared dimension names and
> sizes have nothing to do with the ones use in the grid)
>
> I now get:
>
> $ gdalinfo
> HDF5:"AMSR_U2_L3_SeaIce12km_B04_20120702.he5"://HDFEOS/GRIDS/SpPolarGrid12km/Data_Fields/SI_12km_SH_18H_ASC
> Driver: HDF5Image/HDF5 Dataset
> Files: AMSR_U2_L3_SeaIce12km_B04_20120702.he5
> Size is 632, 664
> Coordinate System is:
> PROJCRS["unnamed",
> BASEGEOGCRS["Unknown datum based upon the custom spheroid",
> DATUM["Not specified (based on custom spheroid)",
> ELLIPSOID["Custom spheroid",6378273,0,
> LENGTHUNIT["metre",1,
> ID["EPSG",9001,
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122,
> CONVERSION["Polar Stereographic (variant B)",
> METHOD["Polar Stereographic (variant B)",
> ID["EPSG",9829]],
> PARAMETER["Latitude of standard parallel",-70,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8832]],
> PARAMETER["Longitude of origin",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8833]],
> PARAMETER["False easting",0,
> LENGTHUNIT["metre",1],
> ID["EPSG",8806]],
> PARAMETER["False northing",0,
> LENGTHUNIT["metre",1],
> ID["EPSG",8807]]],
> CS[Cartesian,2],
> AXIS["(E)",north,
> MERIDIAN[90,
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> ORDER[1],
> LENGTHUNIT["Meter",1]],
> AXIS["(N)",north,
> MERIDIAN[0,
> ANGLEUNIT["degree",0.0174532925199433,
> ID["EPSG",9122]]],
> ORDER[2],
> LENGTHUNIT["Meter",1]]]
> Data axis to CRS axis mapping: 1,2
> Origin = (-395.000,435.000)
> Pixel Size = (12500.000,-12500.000)
> Metadata:
>   Conventions=CF-1.6
>   HDFEOS_INFORMATION_HDFEOSVersion=HDFEOS_5.1.15
>   history=This version of the Sea Ice processing code contains updates
> provided by the science team on September 16, 2019. For details on these
> updates, see the release notes provided in the DAP.
>   institution=NASA's AMSR Science Investigator-led Processing System (SIPS)
>   references=Please cite these data as: Markus, T., J. C. Comiso, and W.
> N. Meier. 2018. AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness
> Temperatures, Sea Ice Concentration, Motion & Snow Depth Polar Grids,
> Version 1. [Indicate subset used]. Boulder, Colorado USA. NASA National
> Snow and Ice Data Center Distributed Active Archive Center. doi:
> https://doi.org/10.5067/RA1MIJOYPK3P.
>   source=satellite observation
>   title=AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures, Sea
> Ice Concentration, Motion & Snow Depth Polar Grids
> Corner Coordinates:
> Upper Left  (-395.000, 435.000) ( 42d14'27.21"W, 39d11'27.54"S)
> Lower Left  (-395.000,-395.000) (135d 0' 0.00"W, 41d23'59.41"S)
> Upper Right ( 395.000, 435.000) ( 42d14'27.21"E, 39d11'27.54"S)
> Lower Right ( 395.000,-395.000) (135d 0' 0.00"E, 41d23'59.41"S)
> Center  (   0.000,  20.000) (  0d 0' 0.01"E, 88d 8'51.76"S

[gdal-dev] HDF5 and geolocation arrays

[Michael Sumner via gdal-dev]

This HDF5 (requires earthdata credentials your "Authorization: Bearer
" in GDAL_HTTP_HEADERS, or equiv) presents without geolocation
arrays.

gdalinfo "/vsicurl/
https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5;
-sd 26
Driver: HDF5Image/HDF5 Dataset
Files: /vsicurl/
https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
Size is 608, 896
Metadata:
  Conventions=CF-1.6
  HDFEOS_INFORMATION_HDFEOSVersion=HDFEOS_5.1.15
  history=This version of the Sea Ice processing code contains updates
provided by the science team on September 16, 2019. For details on these
updates, see the release notes provided in the DAP.
  institution=NASA's AMSR Science Investigator-led Processing System (SIPS)
  references=Please cite these data as: Markus, T., J. C. Comiso, and W. N.
Meier. 2018. AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures,
Sea Ice Concentration, Motion & Snow Depth Polar Grids, Version 1.
[Indicate subset used]. Boulder, Colorado USA. NASA National Snow and Ice
Data Center Distributed Active Archive Center. doi:
https://doi.org/10.5067/RA1MIJOYPK3P.
  source=satellite observation
  title=AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures, Sea
Ice Concentration, Motion & Snow Depth Polar Grids
Corner Coordinates:
Upper Left  (0.0,0.0)
Lower Left  (0.0,  896.0)
Upper Right (  608.0,0.0)
Lower Right (  608.0,  896.0)
Center  (  304.0,  448.0)
Band 1 Block=608x1 Type=Int32, ColorInterp=Undefined
  Metadata:
comment=data value meaning: 0 -- Open Water, 110 -- missing/not
calculated, 120 -- Land
coordinates=lon lat
long_name=Sea ice concentration daily average
units=percent



gdalinfo --version
GDAL 3.9.0dev-cb4d30f56d, released 2024/04/15

The geolocation arrays are sds 33 and 32 respectively:

HDF5:"/vsicurl/
https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
"://HDFEOS/GRIDS/NpPolarGrid12km/lon

HDF5:"/vsicurl/
https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5
"://HDFEOS/GRIDS/NpPolarGrid12km/lat

And things work when lining those up in VRT with warp. Can the HDF5 driver
be made to auto-detect these geolocation arrays?

I see that the NETCDF driver actually does:

gdalinfo "NetCDF:/vsicurl/
https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5;
-sd 26

I'm asking as an email rather than pursuing the fix because, these data are
actually a regular grid on the north and south poles, and so geolocation by
arrays is sub-optimal  the specification is listed in

https://nsidc.org/sites/default/files/au_si12-v001-userguide_1.pdf

and the two parameter sets are

Np-north: -te -385,  -535, 375, 585 -t_srs EPSG:3411
Sp-south: -te -395,  -395, 395, 435 -t_srs EPSG:3412

Is this generally something we should pursue within GDAL? It seems like an
endless task to detect-on-open exactly this situation and assign the easy
fix, but this is a pretty fundamental data stream and it's very common so
the longlat/arrays might be numerically detectable with other
heuristics hinting that it's polar (??) and there are plenty of other
sources that present equivalents in the right way e.g. this one:

"/vsicurl/
https://noaadata.apps.nsidc.org/NOAA/G02135/north/daily/geotiff/2012/07_Jul/N_20120702_concentration_v3.0.tif
"

The right approach is probably to inform the providers and get the right
metadata baked in ... but there's pros and cons to either. I'm not sure
there's even enough information in these files to clearly detect the
situation, it would be a bit like the NSIDCbin driver with its very strict
requirements.

Cheers, Mike


--

Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] lost python api doc?

[Michael Sumner via gdal-dev]

Ah awesome, thank you!



On Mon, Apr 29, 2024 at 8:59 AM Even Rouault 
wrote:

> Michael,
>
> The Python doc got reorganized one month ago (
> https://github.com/OSGeo/gdal/pull/9575), so perhaps Google hasn't
> updated yet the new references
>
> If you look for gdal.warp on gdal.org directly :
> https://gdal.org/search.html?q=gdal.warp_keywords=yes=default
>
> that leads to https://gdal.org/api/python/utilities.html#osgeo.gdal.Warp
>
> Even
> Le 29/04/2024 à 00:52, Michael Sumner via gdal-dev a écrit :
>
> I'm confused about how to browse the python api docs, I used to just web
> search "osgeo.gdal Warp", and then scan down the first result page to find
> "Warp(" and I had the documentation I needed.
>
> Where is that now? (Why can't I find it, sorry - grepping the sources for
> "Warp(" only results in RFC 59 references, so I think this might be a real
> problem).
>
> If I'm just navigating the API docs incorrectly I'd appreciate a
> walkthrough to find the signature and param descripts for "Warp()". Thanks!
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] lost python api doc?

[Michael Sumner via gdal-dev]

I'm confused about how to browse the python api docs, I used to just web
search "osgeo.gdal Warp", and then scan down the first result page to find
"Warp(" and I had the documentation I needed.

Where is that now? (Why can't I find it, sorry - grepping the sources for
"Warp(" only results in RFC 59 references, so I think this might be a real
problem).

If I'm just navigating the API docs incorrectly I'd appreciate a
walkthrough to find the signature and param descripts for "Warp()". Thanks!

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Reading interpolated values on DSM

[Michael Sumner via gdal-dev]

Or a grouping function that returned the cell index for neighbours and
weighting that are involved in whatever calculation summary is wanted.

Maybe the warper could return this as a starting point rather than doing
the "task at hand". ?



On Wed, Apr 24, 2024 at 8:51 PM Even Rouault via gdal-dev <
gdal-dev@lists.osgeo.org> wrote:

> Hi,
>
> I guess this discussion, and past similar ones, calls for an enhancement.
> A new API function, like CPLErr
> GDALRasterInterpolateAtPoint(GDALRasterBandH, double dfPixel, double
> dfLocation, GDALRIOResampleAlg eInterpolation, double *pdfValue), that
> could be used by a new mode of gdallocationinfo. The GDALRPCGetDEMHeight()
> function in alg/gdal_rpc.cpp is a plausible candidate implementation for
> bilinear and bicubic (we could potentially restrict to that at the moment).
>
> Even
> Le 24/04/2024 à 10:33, Javier Jimenez Shaw via gdal-dev a écrit :
>
> Hi
>
> I would like to read in QGIS or GDAL an interpolated value in a DSM (well,
> actually it is a geoid model, but it is the same behaviour). See that I do
> not want the pixel value, but the linear interpolation among the neighbour
> pixels, assuming that the pixel value is in the center of the pixel.
> For instance, this file
> https://www.isgeoid.polimi.it/Geoid/Asia/Japan/japan2000_g.html
>
> Is there any way to get it (without implementing the interpolation myself)?
> If I understood correctly gdallocationinfo is not interpolating, just
> giving the pixel value.
>
> Thanks
>
> .___ ._ ..._ .. . ._.  .___ .. __ . _. . __..  ...  ._ .__
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] advice on python style - get file/dataset as bytes

[Michael Sumner via gdal-dev]

Hi, I'm getting some pushback on my code style. I just want the raw bytes
in-memory of a file in Python, for a manageable tiny dataset.

Is there anything especially wrong with the following?   It's essentially
the same as the answer here:

https://lists.osgeo.org/pipermail/gdal-dev/2016-August/045030.html

   (I know it needs error checks, I just mean the bare style of how to get
the task done).

temp_name = "/vsimem/some.tif"
from osgeo import gdal
gdal.UseExceptions()

## we do something with vrt
dsn = "vrt:///vsicurl/
https://github.com/cran/rgl/raw/master/inst/textures/worldsmall.png?a_ullr=-180,90,180,-90=rgb
"
ds = gdal.Open(dsn)
## write to tif format, but using MEM
gdal.Translate(temp_name, ds)

## now obtain those bytes from the GDAL virtual file
f = gdal.VSIFOpenL(temp_name, 'rb')
gdal.VSIFSeekL(f, 0, 2) # end
size = gdal.VSIFTellL(f)
gdal.VSIFSeekL(f, 0, 0) # begin

## the desired result is this
data = gdal.VSIFReadL(1, size, f)

## and clean up
gdal.VSIFCloseL(f)
gdal.Unlink(temp_name)

Cheers, Mike



-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] discussion on xarray and the missing abstract-grid specification

[Michael Sumner via gdal-dev]

ah thanks Even, I understand that and it was the first example that Ryan
acknowledged was a problem a while back - he'd seen the Float32 problem in
GHRSST. Excellent to get your take and the details.

I've written up a small description of how this affects GHRSST, and my
still-undecided next step to fix it (as COGs).

https://gist.github.com/mdsumner/66aac8e54b54faa24e14b63f396dca32

(this was going to be part of my pitch to xarray and now I don't need so
much extra background to make the point).

I'm keen to collate examples if people have good ones, I have plenty though


Cheers, Mike

On Thu, 4 Apr 2024, 08:26 Even Rouault,  wrote:

> Hi,
>
> Not that I've a strong opinion on what GeoZarr/XArray should do, but yes
> 1D coordinate arrays are a bit of a pain to deal with, when they actually
> just encode a regularly spaced grid. This is something I've hit in the
> bridge between the GDAL multidimensional API (roughly netCDF based) and the
> GDAL "classic" 2D raster API, when you want to expose a view of a multi-dim
> dataset as a classic one. There's this bit of logic in
> https://github.com/OSGeo/gdal/blob/master/gcore/gdalmultidim.cpp#L7334 to
> guess if the 1D coordinate array is regularly spaced or not. The netCDF
> driver has a similar heuristics too. This is quite error prone, especially
> if Float32 coordinates are used. There it can be difficult to detect if it
> is a regularly spaced array but Float32 hasn't sufficient precision to
> fully encode it, or if it there are genuine irregularities in the spacing.
> For Float64, the precision issues are less a practical error, because
> comparing the expected value (assuming a regular spacing) with the one of
> the 1D coordinate array with a very small relative epsilon, like 1e-8, is
> sufficient.
>
> Even
>
>
> Le 03/04/2024 à 22:56, Michael Sumner via gdal-dev a écrit :
>
> Here's an (ahem) extremely important discussion on the prospects for
> xarray to extend from the coordinates-only model (like that of netcdf) for
> geo reference:
>
>
> https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/
>
> I'm heartened to see this recognized at this level. I think GDAL per se
> should feature in this discussion prominently, and not just the downstream
> Python packages that are mentioned.
>
> For my input, I'll be highlighting examples of
>
> -  degenerate rectilinear coordinates, and the problem of whether to
> assign a regular grid there (with a trivial lightweight Translate a la
> -a_ullr), or to define a new regular grid and push it through the Warp api,
> as a dataset that picks up or can be pointed to the 1D coordinate arrays.
> This has been the crux of the issue for me, I get to decisions where it's
> not clear what was intended or what should be done going forward.
>
> - actual curvilinear coordinates, and the very very general power of that
> to resolve to a regular grid with very simple specification (specify some
> or nothing of target/extent/resolution/crs/dimension). Key problems here
> are when the coordinate arrays are not auto-detected (rare) or when the
> longitudes are unwrapped (less rare).
>
> - that extent+dimension is complementary to the affine transform, and way
> more user-friendly way to work with rasters for the assign and/or target
> specification step (this is huge in R since the {raster} package ca. 2009
> and very well supported by Warp and Translate). It's true that
> extent+dimension can't do shear params but it's a perfectly valid way to
> work and flipping from transform to extent is trivial numerically.
>
> I hope some more voices from this community pitch in as well, very happy
> to discuss here or offline about any of this.  I'm not really a technical
> expert but I have a good overview of the landscape that GDAL sits in, and
> I'm making similar pitches for involvement in the R communities.
>
> Cheers, Mike
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 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

[gdal-dev] discussion on xarray and the missing abstract-grid specification

[Michael Sumner via gdal-dev]

Here's an (ahem) extremely important discussion on the prospects for xarray
to extend from the coordinates-only model (like that of netcdf) for geo
reference:

https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/

I'm heartened to see this recognized at this level. I think GDAL per se
should feature in this discussion prominently, and not just the downstream
Python packages that are mentioned.

For my input, I'll be highlighting examples of

-  degenerate rectilinear coordinates, and the problem of whether to assign
a regular grid there (with a trivial lightweight Translate a la -a_ullr),
or to define a new regular grid and push it through the Warp api, as a
dataset that picks up or can be pointed to the 1D coordinate arrays. This
has been the crux of the issue for me, I get to decisions where it's not
clear what was intended or what should be done going forward.

- actual curvilinear coordinates, and the very very general power of that
to resolve to a regular grid with very simple specification (specify some
or nothing of target/extent/resolution/crs/dimension). Key problems here
are when the coordinate arrays are not auto-detected (rare) or when the
longitudes are unwrapped (less rare).

- that extent+dimension is complementary to the affine transform, and way
more user-friendly way to work with rasters for the assign and/or target
specification step (this is huge in R since the {raster} package ca. 2009
and very well supported by Warp and Translate). It's true that
extent+dimension can't do shear params but it's a perfectly valid way to
work and flipping from transform to extent is trivial numerically.

I hope some more voices from this community pitch in as well, very happy to
discuss here or offline about any of this.  I'm not really a technical
expert but I have a good overview of the landscape that GDAL sits in, and
I'm making similar pitches for involvement in the R communities.

Cheers, Mike

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] STACIT and limit

[Michael Sumner via gdal-dev]

Oh, my bad. It's documented to split into SDS for multiple CRS.

Ouch that's not how I thought it was working, but it makes sense.

Thanks, Mike



On Thu, Apr 4, 2024 at 1:26 AM Michael Sumner  wrote:

> this works for me, there are ten items in the filelist:
>
> gdalinfo "STACIT:\"
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z\":asset=visual;
> -oo MAX_ITEMS=10
>
> and with 20 it's also fine, filelist of 20.
>
> but, with MAX_ITEMS=40 there's the warning
>
> Warning 1: Maximum number of items (1000) allowed to be retrieved has been
> hit
>
> and the result is SDS list split by CRS:
>
> gdalinfo "STACIT:\"
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z\":asset=visual;
> -oo MAX_ITEMS=40
> Driver: VRT/Virtual Raster
> Files: none associated
> Size is 0, 0
> Subdatasets:
>   SUBDATASET_1_NAME=STACIT:"
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
> ":asset=visual,crs=EPSG_32631
>   SUBDATASET_1_DESC=Asset visual of
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
> in CRS EPSG:32631
>   SUBDATASET_2_NAME=STACIT:"
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
> ":asset=visual,crs=EPSG_32632
>   SUBDATASET_2_DESC=Asset visual of
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
> in CRS EPSG:32632
>
>
> Another case is this where there should be 4 TCI.tif items, but I only get
> 2 subdatasets (one for each CRS)
>
> gdalinfo "STACIT:\"
> https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=179,-17,180,-16=2023-01-01T00:00:00Z/2023-01-03T23:59:59Z\
> ":asset=visual"
>
> It feels like something awry with the limit for the server, or the driver,
> or maybe interacting.
>
> But, I think it's reasonable expect to not get SDS listed when an asset is
> specified?
>
> Cheers, Mike
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] HDF5 and identified fields / primary dimension

[Michael Sumner via gdal-dev]

> For that particular file, I see that the "feature_id" variable
> (corresponding to the "feature_id" dimension) has a cf_role =
> "timeseries_id" attribute, and that the global metadata has a
> featureType = "timeSeries" attribute. So given
>
> https://cfconventions.org/Data/cf-conventions/cf-conventions-1.11/cf-conventions.html#coordinates-metadata
> , this seems to be relatively standardized, and in that case the
> heuristics could be improve to recognize that the main dimension is
> feature_id (probably with a test that the size of the time dimension is
> 1).  As far as I can see/remember, the vector layer support in netCDF
> was originally developed for the featureType=point and profile use cases
> , so some tuning for timeseries isn't unexpected
>
>
Thanks!  I've made *some* progress, the deepest I've been down in that file
... I hope to be able to craft these suggestions at some point.

Cheers, Mike



> Even
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] STACIT and limit

[Michael Sumner via gdal-dev]

this works for me, there are ten items in the filelist:

gdalinfo "STACIT:\"
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z\":asset=visual;
-oo MAX_ITEMS=10

and with 20 it's also fine, filelist of 20.

but, with MAX_ITEMS=40 there's the warning

Warning 1: Maximum number of items (1000) allowed to be retrieved has been
hit

and the result is SDS list split by CRS:

gdalinfo "STACIT:\"
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z\":asset=visual;
-oo MAX_ITEMS=40
Driver: VRT/Virtual Raster
Files: none associated
Size is 0, 0
Subdatasets:
  SUBDATASET_1_NAME=STACIT:"
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
":asset=visual,crs=EPSG_32631
  SUBDATASET_1_DESC=Asset visual of
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
in CRS EPSG:32631
  SUBDATASET_2_NAME=STACIT:"
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
":asset=visual,crs=EPSG_32632
  SUBDATASET_2_DESC=Asset visual of
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=0,0,10,10=2023-01-01T00:00:00Z/2023-12-31T23:59:59Z
in CRS EPSG:32632


Another case is this where there should be 4 TCI.tif items, but I only get
2 subdatasets (one for each CRS)

gdalinfo "STACIT:\"
https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a=179,-17,180,-16=2023-01-01T00:00:00Z/2023-01-03T23:59:59Z\
":asset=visual"

It feels like something awry with the limit for the server, or the driver,
or maybe interacting.

But, I think it's reasonable expect to not get SDS listed when an asset is
specified?

Cheers, Mike

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] HDF5 and identified fields / primary dimension

[Michael Sumner via gdal-dev]

well actually, I think what I'm asking for is the intended behaviour, but
there's an error.

Is it meant to detect sets of variables on 1D dimensions and present them
as layers? That's what would make sense to me.

Still exploring.

Cheers, Mike



On Tue, Apr 2, 2024 at 5:36 AM Michael Sumner  wrote:

> This source has an array on 'feature_id' with 2729077 values, with various
> fields
>
> elevation, longitude, latitude, qBtmVertRunoff, qBucket, etc
>
>
> '/vsis3/noaa-nwm-retro-v2.0-pds/full_physics/2017/20170401.CHRTOUT_DOMAIN1.comp'
>
> It is accessible via the mdim api.
>
> Structurally it is basically a table with rows per feature_id and  columns
> per fields, but it has a length-1 pair of fields "time" and
> "reference_time" defined on dimension time, this is like a single time step
> per file (like an unlimited dimension in the classic 2D case).
>
> Accessing with the vector API reports that it can't treat this as a table
> because of those time values that don't match the feature_id dimension:
>
> ogrinfo
> NETCDF:'/vsis3/noaa-nwm-retro-v2.0-pds/full_physics/2017/20170401.CHRTOUT_DOMAIN1.comp'
> -ro
>
> Warning 1: The dataset has several variables that could be identified as
> vector fields, but not all share the same primary dimension. Consequently
> they will be ignored.
>
> I've seen similar cases in other files. I presume the driver could be
> updated to 1) choose the primary dimension and read the values while ignore
> others 2) user-specify the dimension to include, or 3) user-specify the
> fields to exclude
>
> So:
>
> - is there a workaround to enable the vector driver to focus on the
> primary dimension?
> - would a PR along those lines have to consider greater difficulties than
> applying the proposed updates to arrays using the primary dimension only?
>  I'd only consider this for strictly 1D arrays.
> - degenerate dimensions could be used to copy-out the value of the other
> dims (I'd consider this an optional extra)
>
> (It's a bit special-case-y, you wouldn't want to go to multi-arrays and
> have them flatten out multi-dims in a general way, I think, but degenerate
> dimensions might be worth consideration )
>
> Appreciate any thoughts, thanks! I'd quite like to have the
> vector-approach work as well as the mdim approach, I think they are nicely
> complementary and provide different pros and cons.
>
> Cheers, Mike
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] HDF5 and identified fields / primary dimension

[Michael Sumner via gdal-dev]

This source has an array on 'feature_id' with 2729077 values, with various
fields

elevation, longitude, latitude, qBtmVertRunoff, qBucket, etc

'/vsis3/noaa-nwm-retro-v2.0-pds/full_physics/2017/20170401.CHRTOUT_DOMAIN1.comp'

It is accessible via the mdim api.

Structurally it is basically a table with rows per feature_id and  columns
per fields, but it has a length-1 pair of fields "time" and
"reference_time" defined on dimension time, this is like a single time step
per file (like an unlimited dimension in the classic 2D case).

Accessing with the vector API reports that it can't treat this as a table
because of those time values that don't match the feature_id dimension:

ogrinfo
NETCDF:'/vsis3/noaa-nwm-retro-v2.0-pds/full_physics/2017/20170401.CHRTOUT_DOMAIN1.comp'
-ro

Warning 1: The dataset has several variables that could be identified as
vector fields, but not all share the same primary dimension. Consequently
they will be ignored.

I've seen similar cases in other files. I presume the driver could be
updated to 1) choose the primary dimension and read the values while ignore
others 2) user-specify the dimension to include, or 3) user-specify the
fields to exclude

So:

- is there a workaround to enable the vector driver to focus on the primary
dimension?
- would a PR along those lines have to consider greater difficulties than
applying the proposed updates to arrays using the primary dimension only?
 I'd only consider this for strictly 1D arrays.
- degenerate dimensions could be used to copy-out the value of the other
dims (I'd consider this an optional extra)

(It's a bit special-case-y, you wouldn't want to go to multi-arrays and
have them flatten out multi-dims in a general way, I think, but degenerate
dimensions might be worth consideration )

Appreciate any thoughts, thanks! I'd quite like to have the vector-approach
work as well as the mdim approach, I think they are nicely complementary
and provide different pros and cons.

Cheers, Mike

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] palette file format doc ?

[Michael Sumner via gdal-dev]

>
>
>
> > And, can index be *value* in any contexts?
>
> If you use a raster with a signed data type, that could be negative
> values (assuming I understand your question)
>
>

ah I see, arbitrary integer values map to a colour 0:(n-1) colours, match
the ordered n values in the raster - that is "value" in the way I was
thinking  ( I had in mind more general cases like interval ranges, but
that's  another level up)

I might be entirely messing up the concept ... will check my take by
experiment.

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] palette file format doc ?

[Michael Sumner via gdal-dev]

Is the palette_file  .txt format documented?

https://gdal.org/programs/gdalattachpct.html

It's mentioned in a few utilities, and created by tests but I couldn't find
an existing example or a description (I guessed, incorrectly at first,
leaving out the index column). I take it that it is 0-255 range values, but
is index 0-based or 1? And, can index be *value* in any contexts?

index  red green blue 

If it's documented I'll link it in PRs. I see qml and qlr here:
https://docs.qgis.org/3.34/en/docs/user_manual/appendices/qgis_file_formats.html

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] overview sizes for COG

[Michael Sumner via gdal-dev]

Excellent!  Thanks for the answers.

When I've explored a bit more I might implement the overview sizes just so
we can match downstream tools (the current motivation is like-for-like
performance comparison, I haven't looked there but I think odc goes very
low level to eke out speed).

Cheers, Mike



On Wed, Mar 20, 2024 at 10:39 AM Even Rouault 
wrote:

> Michael,
>
> Hi, can we specify overview sizes exactly?
>
> No, the BuildOverviews() interface onlys take an array of overview
> factors. It is up to the driver implementation to decide how it computes
> the overview size from the main raster size and overview factor.
>
> The COG driver is a bit of a special case, since it automagically compute
> overviews if they don't exist in the source dataset, and doesn't take an
> explicit list of overview factors. It simply divides by two successively,
> with floor truncation:
> https://github.com/OSGeo/gdal/blob/master/frmts/gtiff/cogdriver.cpp#L1092
>
> I have this odd grid that is 36000x17999, and I get consequently yucky
> overview sizes:
>
> gdal_create -outsize 36000 17999 -ot Int8 -co SPARSE_OK=YES -a_srs
> EPSG:4326 -a_ullr 0 17999 0 36000 weird.tif
> gdal_translate weird.tif cog.tif  -of COG
> gdalinfo cog.tif | grep Overviews
> #  Overviews: 18000x8999, 9000x4499, 4500x2249, 2250x1124, 1125x562,
> 562x281, 281x14
>
> Should I worry about this?
>
> Probably not ? :-)
>
> Can the sizes be specified via Translate? (I can see how they could be
> generated externally and bundled, but I'm looking for a way at this high
> level).
>
> I guess someone could potentially add a OVERVIEW_SIZES=W1xH1,W2xH2,...
> creation option to the COG driver if that was really needed...
>
>
> ( Maybe I shouldn't even care about the sizes ... and maybe I should
> resample the grid to be a better overall size, but I consider that out of
> scope here.)
>
> I'm comparing with other COG creators, e.g. odc which makes tidy overviews
> ( I will explore how it's getting done).
>
> It looks like they might use (previous_size + 1) / 2 (or ceil rounding)
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] overview sizes for COG

[Michael Sumner via gdal-dev]

Hi, can we specify overview sizes exactly? I have this odd grid that is
36000x17999, and I get consequently yucky overview sizes:

gdal_create -outsize 36000 17999 -ot Int8 -co SPARSE_OK=YES -a_srs
EPSG:4326 -a_ullr 0 17999 0 36000 weird.tif
gdal_translate weird.tif cog.tif  -of COG
gdalinfo cog.tif | grep Overviews
#  Overviews: 18000x8999, 9000x4499, 4500x2249, 2250x1124, 1125x562,
562x281, 281x14

Should I worry about this?  Can the sizes be specified via Translate? (I
can see how they could be generated externally and bundled, but I'm looking
for a way at this high level).

( Maybe I shouldn't even care about the sizes ... and maybe I should
resample the grid to be a better overall size, but I consider that out of
scope here.)

I'm comparing with other COG creators, e.g. odc which makes tidy overviews
( I will explore how it's getting done).

import xarray
from odc.geo.xr import assign_crs
from odc.geo.cog import write_cog

ds = xarray.open_dataset("weird.tif")
ds = assign_crs(ds, crs="EPSG:4326")  ## do something to tag it as odc


write_cog(ds["band_data"], "odc.tif")
os.system("gdalinfo odc.tif | grep Overviews")
#  Overviews: 18000x9000, 9000x4500, 4500x2250, 2250x1125, 1125x563

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] netcdf url : should lack of range support be the error?

[Michael Sumner via gdal-dev]

on this URL I get an error, and as far as i understand - the server doesn't
support range downloading, the file otherwise works fine.

gdalinfo failed - unable to open '/vsicurl/
https://erddap.emodnet.eu/erddap/files/biology_6640_benthos_NorthSea_e4af_0f0e_6a73/04_2021_6640_diva_benthos_erddap.nc
'.


Should this trigger "Range downloading is not supported by this server" ?

Full output with debug messages on:


gdalinfo --config CPL_DEBUG ON /vsicurl/
https://erddap.emodnet.eu/erddap/files/biology_6640_benthos_NorthSea_e4af_0f0e_6a73/04_2021_6640_diva_benthos_erddap.nc
-if NetCDF


HTTP: libcurl/7.81.0 OpenSSL/3.0.2 zlib/1.2.11 brotli/1.0.9 zstd/1.4.8
libidn2/2.3.2 libpsl/0.21.0 (+libidn2/2.3.2) libssh/0.9.6/openssl/zlib
nghttp2/1.43.0 librtmp/2.3 OpenLDAP/2.5.16
VSICURL: GetFileSize(
https://erddap.emodnet.eu/erddap/files/biology_6640_benthos_NorthSea_e4af_0f0e_6a73/04_2021_6640_diva_benthos_erddap.nc)=339890640
 response_code=200
VSICURL: Downloading 0-16383 (
https://erddap.emodnet.eu/erddap/files/biology_6640_benthos_NorthSea_e4af_0f0e_6a73/04_2021_6640_diva_benthos_erddap.nc).
..
VSICURL: Got response_code=416
ERROR 4: `/vsicurl/
https://erddap.emodnet.eu/erddap/files/biology_6640_benthos_NorthSea_e4af_0f0e_6a73/04_2021_6640_diva_benthos_erddap.nc'
not recognized as being in a supported file format.


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

I don't understand how jammy is "old" when the full build is itself using
"BASE_IMAGE=ubuntu:22.04"

But, I'm out of my depth in these emails and trying to learn, thanks!

Cheers, Mike




On Wed, 7 Feb 2024, 00:46 Javier Jimenez Shaw,  wrote:

> Could you set up your VMs to include those SSE instructions? I think that
> keeping VMs that "old" configured is a source of problems using
> pre-compiled binaries.
> The same way GDAL updates dependencies of compilers and other libraries to
> something more modern (but not too modern), those SSE instructions should
> be updated.
>
> @Even knowing now that the "old" hardware is "virtually old", should we
> remove AVX2 compatibility from ubuntu-full-latest? I do not know how much
> is the performance impact.
>
> Cheers,
> Javier.
>
> On Mon, 5 Feb 2024 at 21:54, Michael Sumner  wrote:
>
>> yes, jammy VM on openstack is the host (and is where I run pretty much
>> everything, though will increasingly use AWS).
>>
>> Thanks for the note, I'll try on other systems too. We need a
>> security-allow set for vsicurl to work so if there are other little details
>> I'll be keen to flush them out.
>>
>> Cheers, Mike
>>
>> On Mon, 5 Feb 2024, 18:44 Javier Jimenez Shaw, 
>> wrote:
>>
>>> Hi Mike
>>>
>>> Out of curiosity, are run running it in a virtual machine?
>>> A few year ago I had problems running a program in a virtual machine
>>> (virtualbox, but I read it happens in others) due to a missing SSE
>>> instruction. The solution there was to "enable" the missing instructions in
>>> the virtual machine configuration (that I don't know why it was not the
>>> default).
>>>
>>> Cheers
>>>
>>>
>>> On Mon, 5 Feb 2024, 02:04 Even Rouault via gdal-dev, <
>>> gdal-dev@lists.osgeo.org> wrote:
>>>
>>>> ghcr.io/osgeo/gdal:ubuntu-full-latest has been regenerated with the
>>>> rebuild of TileDB without AVX2. I've also enabled the
>>>> drivers-with-external-depencies-built-as-plugin GDAL build mode, so it is
>>>> easy to just remove a given plugin by deleting the corresponding .so in
>>>> /usr/lib/x86_64-linux-gnu/gdalplugins
>>>>
>>>> Even
>>>> Le 04/02/2024 à 22:51, Michael Sumner a écrit :
>>>>
>>>> indeed there's no avx2:
>>>>
>>>> cat /proc/cpuinfo|grep sse|head -n 1
>>>> flags   : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge
>>>> mca cmov pat pse36 clflush mmx fxsr sse sse2 syscall nx mmxext fxsr_opt
>>>> pdpe1gb rdtscp lm rep_good nopl cpuid extd_apicid tsc_known_freq pni
>>>> pclmulqdq ssse3 fma cx16 sse4_1 sse4_2 x2apic popcnt tsc_deadline_timer aes
>>>> xsave avx f16c hypervisor lahf_lm cmp_legacy svm cr8_legacy abm sse4a
>>>> misalignsse 3dnowprefetch osvw xop fma4 tbm perfctr_core ssbd ibpb vmmcall
>>>> tsc_adjust bmi1 virt_ssbd arat npt nrip_save arch_capabilities
>>>>
>>>> Cheers, Mike
>>>>
>>>>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

yes, jammy VM on openstack is the host (and is where I run pretty much
everything, though will increasingly use AWS).

Thanks for the note, I'll try on other systems too. We need a
security-allow set for vsicurl to work so if there are other little details
I'll be keen to flush them out.

Cheers, Mike

On Mon, 5 Feb 2024, 18:44 Javier Jimenez Shaw,  wrote:

> Hi Mike
>
> Out of curiosity, are run running it in a virtual machine?
> A few year ago I had problems running a program in a virtual machine
> (virtualbox, but I read it happens in others) due to a missing SSE
> instruction. The solution there was to "enable" the missing instructions in
> the virtual machine configuration (that I don't know why it was not the
> default).
>
> Cheers
>
>
> On Mon, 5 Feb 2024, 02:04 Even Rouault via gdal-dev, <
> gdal-dev@lists.osgeo.org> wrote:
>
>> ghcr.io/osgeo/gdal:ubuntu-full-latest has been regenerated with the
>> rebuild of TileDB without AVX2. I've also enabled the
>> drivers-with-external-depencies-built-as-plugin GDAL build mode, so it is
>> easy to just remove a given plugin by deleting the corresponding .so in
>> /usr/lib/x86_64-linux-gnu/gdalplugins
>>
>> Even
>> Le 04/02/2024 à 22:51, Michael Sumner a écrit :
>>
>> indeed there's no avx2:
>>
>> cat /proc/cpuinfo|grep sse|head -n 1
>> flags   : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge
>> mca cmov pat pse36 clflush mmx fxsr sse sse2 syscall nx mmxext fxsr_opt
>> pdpe1gb rdtscp lm rep_good nopl cpuid extd_apicid tsc_known_freq pni
>> pclmulqdq ssse3 fma cx16 sse4_1 sse4_2 x2apic popcnt tsc_deadline_timer aes
>> xsave avx f16c hypervisor lahf_lm cmp_legacy svm cr8_legacy abm sse4a
>> misalignsse 3dnowprefetch osvw xop fma4 tbm perfctr_core ssbd ibpb vmmcall
>> tsc_adjust bmi1 virt_ssbd arat npt nrip_save arch_capabilities
>>
>> Cheers, Mike
>>
>>
>>
>> On Sun, Feb 4, 2024 at 10:55 PM Even Rouault 
>> wrote:
>>
>>> ok, so I believe this is the AVX2 issue I was talking about, as I
>>> realize that enabling AVX2 is the default mode when TileDB is built from
>>> source (which the Docker image does), and must be explicitly disabled with
>>> "./bootstrap --disable-avx2" (I've just changed the build recipe to include
>>> that, will take effect next time the images are refreshed)
>>>
>>> To confirm, can you send or just check the output of : cat
>>> /proc/cpuinfo|grep sse|head -n 1
>>>
>>> If there is no "avx2" in it, this is at 99.9% the reason of the issue.
>>>
>>> Even
>>> Le 04/02/2024 à 06:20, Michael Sumner a écrit :
>>>
>>> skipping TileDB does fix:
>>>
>>> ogr2ogr /tmp/newdir
>>> https://github.com/SymbolixAU/geojsonsf/raw/master/inst/examples/geo_melbourne.geojson
>>>  -f
>>> "ESRI Shapefile"
>>> export GDAL_SKIP=TileDB
>>> ogrinfo /tmp/newdir/
>>> INFO: Open of `/tmp/newdir/'
>>>   using driver `ESRI Shapefile' successful.
>>> 1: geo_melbourne (Polygon)
>>>
>>> unset GDAL_SKIP
>>> ogrinfo /tmp/newdir/
>>> Illegal instruction (core dumped)
>>>
>>> I failed to explain that I'm using gdal containers from the repo:
>>>
>>> docker run --rm -ti ghcr.io/osgeo/gdal:ubunt
>>>
>>> u-full-latest
>>>
>>> apt update
>>> apt install -y gdb
>>>
>>> Here's the output of under gdb as you suggested, there was a lot so I
>>> put it on a gist:
>>> https://gist.github.com/mdsumner/839ae6e05ededf640b65bfee3a20a4c0
>>>
>>> gdb --args ogrinfo /tmp/newdir/
>>> > run
>>> > thread apply all bt
>>>
>>> Thanks!
>>>
>>>
>>>
>>>
>>>
>>> On Sat, Feb 3, 2024 at 7:49 PM Even Rouault 
>>> wrote:
>>>
>>>> - When it crashes under gdb, type "thread apply all bt" to get the
>>>> stack trace of all threads
>>>>
>>>> - I suspect there is a connection with
>>>> https://github.com/OSGeo/gdal/pull/9170 , but that pull request
>>>> wouldn't help here as "/tmp/newdir" could be a valid connection to TileDB
>>>>
>>>> - how did you get TileDB installed? It looks to be packaged? Which
>>>> distribution do you use?
>>>>
>>>> - SIGILL reminds me of issues with some TileDB builds using the AVX2
>>>> instruction set by default, which could cause some crash on host CPUs that
&g

Re: [gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

Excellent, thanks Even - confirmed at my end and on my original use-case.

Cheers, Mike



On Mon, Feb 5, 2024 at 12:04 PM Even Rouault 
wrote:

> ghcr.io/osgeo/gdal:ubuntu-full-latest has been regenerated with the
> rebuild of TileDB without AVX2. I've also enabled the
> drivers-with-external-depencies-built-as-plugin GDAL build mode, so it is
> easy to just remove a given plugin by deleting the corresponding .so in
> /usr/lib/x86_64-linux-gnu/gdalplugins
>
> Even
> Le 04/02/2024 à 22:51, Michael Sumner a écrit :
>
> indeed there's no avx2:
>
> cat /proc/cpuinfo|grep sse|head -n 1
> flags   : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca
> cmov pat pse36 clflush mmx fxsr sse sse2 syscall nx mmxext fxsr_opt pdpe1gb
> rdtscp lm rep_good nopl cpuid extd_apicid tsc_known_freq pni pclmulqdq
> ssse3 fma cx16 sse4_1 sse4_2 x2apic popcnt tsc_deadline_timer aes xsave avx
> f16c hypervisor lahf_lm cmp_legacy svm cr8_legacy abm sse4a misalignsse
> 3dnowprefetch osvw xop fma4 tbm perfctr_core ssbd ibpb vmmcall tsc_adjust
> bmi1 virt_ssbd arat npt nrip_save arch_capabilities
>
> Cheers, Mike
>
>
>
> On Sun, Feb 4, 2024 at 10:55 PM Even Rouault 
> wrote:
>
>> ok, so I believe this is the AVX2 issue I was talking about, as I realize
>> that enabling AVX2 is the default mode when TileDB is built from source
>> (which the Docker image does), and must be explicitly disabled with
>> "./bootstrap --disable-avx2" (I've just changed the build recipe to include
>> that, will take effect next time the images are refreshed)
>>
>> To confirm, can you send or just check the output of : cat
>> /proc/cpuinfo|grep sse|head -n 1
>>
>> If there is no "avx2" in it, this is at 99.9% the reason of the issue.
>>
>> Even
>> Le 04/02/2024 à 06:20, Michael Sumner a écrit :
>>
>> skipping TileDB does fix:
>>
>> ogr2ogr /tmp/newdir
>> https://github.com/SymbolixAU/geojsonsf/raw/master/inst/examples/geo_melbourne.geojson
>>  -f
>> "ESRI Shapefile"
>> export GDAL_SKIP=TileDB
>> ogrinfo /tmp/newdir/
>> INFO: Open of `/tmp/newdir/'
>>   using driver `ESRI Shapefile' successful.
>> 1: geo_melbourne (Polygon)
>>
>> unset GDAL_SKIP
>> ogrinfo /tmp/newdir/
>> Illegal instruction (core dumped)
>>
>> I failed to explain that I'm using gdal containers from the repo:
>>
>> docker run --rm -ti ghcr.io/osgeo/gdal:ubunt
>>
>> u-full-latest
>>
>> apt update
>> apt install -y gdb
>>
>> Here's the output of under gdb as you suggested, there was a lot so I put
>> it on a gist:
>> https://gist.github.com/mdsumner/839ae6e05ededf640b65bfee3a20a4c0
>>
>> gdb --args ogrinfo /tmp/newdir/
>> > run
>> > thread apply all bt
>>
>> Thanks!
>>
>>
>>
>>
>>
>> On Sat, Feb 3, 2024 at 7:49 PM Even Rouault 
>> wrote:
>>
>>> - When it crashes under gdb, type "thread apply all bt" to get the stack
>>> trace of all threads
>>>
>>> - I suspect there is a connection with
>>> https://github.com/OSGeo/gdal/pull/9170 , but that pull request
>>> wouldn't help here as "/tmp/newdir" could be a valid connection to TileDB
>>>
>>> - how did you get TileDB installed? It looks to be packaged? Which
>>> distribution do you use?
>>>
>>> - SIGILL reminds me of issues with some TileDB builds using the AVX2
>>> instruction set by default, which could cause some crash on host CPUs that
>>> don't have AVX2 (unlikely on recent hardware though)
>>>
>>> - Setting GDAL_SKIP=TileDB should be a workaround
>>>
>>>
>>> Le 03/02/2024 à 07:15, Michael Sumner a écrit :
>>>
>>> Thanks Even, so there's something about tiledb under gdb (or maybe I am
>>> mangling the context,  I will try variants of the host I'm using). Run with
>>> valgrind included below.
>>>
>>> gdb --args ogrinfo /tmp/newdir/
>>> ...
>>> (gdb) run
>>> Starting program: /usr/local/bin/ogrinfo /tmp/newdir/
>>> [Thread debugging using libthread_db enabled]
>>> Using host libthread_db library
>>> "/lib/x86_64-linux-gnu/libthread_db.so.1".
>>> [New Thread 0x7fffe7757640 (LWP 988)]
>>> [New Thread 0x7fffe6f56640 (LWP 989)]
>>> [New Thread 0x7fffde755640 (LWP 990)]
>>> [New Thread 0x7fffd5f54640 (LWP 991)]
>>> [New Thread 0x7fffc5753640 (LWP 992)]
>>> [New Thread 0x7fffc4f52640 (LWP 993)]
>>> [New Thread 0x7fffb4751640 (LWP 9

Re: [gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

indeed there's no avx2:

cat /proc/cpuinfo|grep sse|head -n 1
flags   : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca
cmov pat pse36 clflush mmx fxsr sse sse2 syscall nx mmxext fxsr_opt pdpe1gb
rdtscp lm rep_good nopl cpuid extd_apicid tsc_known_freq pni pclmulqdq
ssse3 fma cx16 sse4_1 sse4_2 x2apic popcnt tsc_deadline_timer aes xsave avx
f16c hypervisor lahf_lm cmp_legacy svm cr8_legacy abm sse4a misalignsse
3dnowprefetch osvw xop fma4 tbm perfctr_core ssbd ibpb vmmcall tsc_adjust
bmi1 virt_ssbd arat npt nrip_save arch_capabilities

Cheers, Mike



On Sun, Feb 4, 2024 at 10:55 PM Even Rouault 
wrote:

> ok, so I believe this is the AVX2 issue I was talking about, as I realize
> that enabling AVX2 is the default mode when TileDB is built from source
> (which the Docker image does), and must be explicitly disabled with
> "./bootstrap --disable-avx2" (I've just changed the build recipe to include
> that, will take effect next time the images are refreshed)
>
> To confirm, can you send or just check the output of : cat
> /proc/cpuinfo|grep sse|head -n 1
>
> If there is no "avx2" in it, this is at 99.9% the reason of the issue.
>
> Even
> Le 04/02/2024 à 06:20, Michael Sumner a écrit :
>
> skipping TileDB does fix:
>
> ogr2ogr /tmp/newdir
> https://github.com/SymbolixAU/geojsonsf/raw/master/inst/examples/geo_melbourne.geojson
>  -f
> "ESRI Shapefile"
> export GDAL_SKIP=TileDB
> ogrinfo /tmp/newdir/
> INFO: Open of `/tmp/newdir/'
>   using driver `ESRI Shapefile' successful.
> 1: geo_melbourne (Polygon)
>
> unset GDAL_SKIP
> ogrinfo /tmp/newdir/
> Illegal instruction (core dumped)
>
> I failed to explain that I'm using gdal containers from the repo:
>
> docker run --rm -ti ghcr.io/osgeo/gdal:ubunt
>
> u-full-latest
>
> apt update
> apt install -y gdb
>
> Here's the output of under gdb as you suggested, there was a lot so I put
> it on a gist:
> https://gist.github.com/mdsumner/839ae6e05ededf640b65bfee3a20a4c0
>
> gdb --args ogrinfo /tmp/newdir/
> > run
> > thread apply all bt
>
> Thanks!
>
>
>
>
>
> On Sat, Feb 3, 2024 at 7:49 PM Even Rouault 
> wrote:
>
>> - When it crashes under gdb, type "thread apply all bt" to get the stack
>> trace of all threads
>>
>> - I suspect there is a connection with
>> https://github.com/OSGeo/gdal/pull/9170 , but that pull request wouldn't
>> help here as "/tmp/newdir" could be a valid connection to TileDB
>>
>> - how did you get TileDB installed? It looks to be packaged? Which
>> distribution do you use?
>>
>> - SIGILL reminds me of issues with some TileDB builds using the AVX2
>> instruction set by default, which could cause some crash on host CPUs that
>> don't have AVX2 (unlikely on recent hardware though)
>>
>> - Setting GDAL_SKIP=TileDB should be a workaround
>>
>>
>> Le 03/02/2024 à 07:15, Michael Sumner a écrit :
>>
>> Thanks Even, so there's something about tiledb under gdb (or maybe I am
>> mangling the context,  I will try variants of the host I'm using). Run with
>> valgrind included below.
>>
>> gdb --args ogrinfo /tmp/newdir/
>> ...
>> (gdb) run
>> Starting program: /usr/local/bin/ogrinfo /tmp/newdir/
>> [Thread debugging using libthread_db enabled]
>> Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
>> [New Thread 0x7fffe7757640 (LWP 988)]
>> [New Thread 0x7fffe6f56640 (LWP 989)]
>> [New Thread 0x7fffde755640 (LWP 990)]
>> [New Thread 0x7fffd5f54640 (LWP 991)]
>> [New Thread 0x7fffc5753640 (LWP 992)]
>> [New Thread 0x7fffc4f52640 (LWP 993)]
>> [New Thread 0x7fffb4751640 (LWP 994)]
>> [New Thread 0x7fffabf50640 (LWP 995)]
>> [New Thread 0x7fffab74f640 (LWP 996)]
>> [New Thread 0x7fffa2f4e640 (LWP 997)]
>> [New Thread 0x7fff9a74d640 (LWP 998)]
>> [New Thread 0x7fff91f4c640 (LWP 999)]
>> [New Thread 0x7fff8974b640 (LWP 1000)]
>> [New Thread 0x7fff78f4a640 (LWP 1001)]
>> [New Thread 0x7fff78749640 (LWP 1002)]
>> [New Thread 0x7fff6f5ff640 (LWP 1003)]
>>
>> Thread 1 "ogrinfo" received signal SIGILL, Illegal instruction.
>> 0x73773c9e in tiledb::common::ThreadPool::ThreadPool(unsigned
>> long) () from /lib/x86_64-linux-gnu/libtiledb.so.2.16
>>
>>
>>
>>
>> valgrind -s ogrinfo /tmp/newdir
>> ==704== Memcheck, a memory error detector
>> ==704== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
>> ==704== Using Valgrind-3.18.1 and LibVEX; rerun with -h for copyright info
>> ==704== Command: ogrinfo /tmp/newdir
>> ==70

Re: [gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

skipping TileDB does fix:

ogr2ogr /tmp/newdir
https://github.com/SymbolixAU/geojsonsf/raw/master/inst/examples/geo_melbourne.geojson
-f
"ESRI Shapefile"
export GDAL_SKIP=TileDB
ogrinfo /tmp/newdir/
INFO: Open of `/tmp/newdir/'
  using driver `ESRI Shapefile' successful.
1: geo_melbourne (Polygon)

unset GDAL_SKIP
ogrinfo /tmp/newdir/
Illegal instruction (core dumped)

I failed to explain that I'm using gdal containers from the repo:

docker run --rm -ti ghcr.io/osgeo/gdal:ubunt

u-full-latest

apt update
apt install -y gdb

Here's the output of under gdb as you suggested, there was a lot so I put
it on a gist:
https://gist.github.com/mdsumner/839ae6e05ededf640b65bfee3a20a4c0

gdb --args ogrinfo /tmp/newdir/
> run
> thread apply all bt

Thanks!





On Sat, Feb 3, 2024 at 7:49 PM Even Rouault 
wrote:

> - When it crashes under gdb, type "thread apply all bt" to get the stack
> trace of all threads
>
> - I suspect there is a connection with
> https://github.com/OSGeo/gdal/pull/9170 , but that pull request wouldn't
> help here as "/tmp/newdir" could be a valid connection to TileDB
>
> - how did you get TileDB installed? It looks to be packaged? Which
> distribution do you use?
>
> - SIGILL reminds me of issues with some TileDB builds using the AVX2
> instruction set by default, which could cause some crash on host CPUs that
> don't have AVX2 (unlikely on recent hardware though)
>
> - Setting GDAL_SKIP=TileDB should be a workaround
>
>
> Le 03/02/2024 à 07:15, Michael Sumner a écrit :
>
> Thanks Even, so there's something about tiledb under gdb (or maybe I am
> mangling the context,  I will try variants of the host I'm using). Run with
> valgrind included below.
>
> gdb --args ogrinfo /tmp/newdir/
> ...
> (gdb) run
> Starting program: /usr/local/bin/ogrinfo /tmp/newdir/
> [Thread debugging using libthread_db enabled]
> Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
> [New Thread 0x7fffe7757640 (LWP 988)]
> [New Thread 0x7fffe6f56640 (LWP 989)]
> [New Thread 0x7fffde755640 (LWP 990)]
> [New Thread 0x7fffd5f54640 (LWP 991)]
> [New Thread 0x7fffc5753640 (LWP 992)]
> [New Thread 0x7fffc4f52640 (LWP 993)]
> [New Thread 0x7fffb4751640 (LWP 994)]
> [New Thread 0x7fffabf50640 (LWP 995)]
> [New Thread 0x7fffab74f640 (LWP 996)]
> [New Thread 0x7fffa2f4e640 (LWP 997)]
> [New Thread 0x7fff9a74d640 (LWP 998)]
> [New Thread 0x7fff91f4c640 (LWP 999)]
> [New Thread 0x7fff8974b640 (LWP 1000)]
> [New Thread 0x7fff78f4a640 (LWP 1001)]
> [New Thread 0x7fff78749640 (LWP 1002)]
> [New Thread 0x7fff6f5ff640 (LWP 1003)]
>
> Thread 1 "ogrinfo" received signal SIGILL, Illegal instruction.
> 0x73773c9e in tiledb::common::ThreadPool::ThreadPool(unsigned
> long) () from /lib/x86_64-linux-gnu/libtiledb.so.2.16
>
>
>
>
> valgrind -s ogrinfo /tmp/newdir
> ==704== Memcheck, a memory error detector
> ==704== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
> ==704== Using Valgrind-3.18.1 and LibVEX; rerun with -h for copyright info
> ==704== Command: ogrinfo /tmp/newdir
> ==704==
> INFO: Open of `/tmp/newdir'
>   using driver `ESRI Shapefile' successful.
> 1: geo_melbourne (Polygon)
> ==704==
> ==704== HEAP SUMMARY:
> ==704== in use at exit: 25,486 bytes in 216 blocks
> ==704==   total heap usage: 15,761 allocs, 15,545 frees, 2,390,169 bytes
> allocated
> ==704==
> ==704== LEAK SUMMARY:
> ==704==definitely lost: 0 bytes in 0 blocks
> ==704==indirectly lost: 0 bytes in 0 blocks
> ==704==  possibly lost: 544 bytes in 1 blocks
> ==704==still reachable: 24,942 bytes in 215 blocks
> ==704== suppressed: 0 bytes in 0 blocks
> ==704== Rerun with --leak-check=full to see details of leaked memory
> ==704==
> ==704== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)
>
>
> ogrinfo /tmp/newdir
> Illegal instruction (core dumped)
>
> Cheers, Mike
>
>
>
>
> On Sat, Feb 3, 2024 at 12:46 PM Even Rouault 
> wrote:
>
>> Michael,
>>
>> I'm wondering if there not might be something wrong with your build or
>> runtime environment. Or there's something subtle, because that works fine
>> for me with my dev build or in the ghcr.io/osgeo/gdal:alpine-normal-3.8.3
>> Docker image
>>
>> Try running "valgrind ogrinfo /tmp/newdir/" or "gdb --args ogrinfo
>> /tmp/newdir/" (type "run") to get more useful information
>>
>> Even
>> Le 03/02/2024 à 02:35, Michael Sumner via gdal-dev a écrit :
>>
>> I'm getting Illegal instruction / core dumped on ogrinfo of a directory:
>>
>> ogr2ogr /tmp/newdir
>> https://github.com/SymbolixAU/g

Re: [gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

Thanks Even, so there's something about tiledb under gdb (or maybe I am
mangling the context,  I will try variants of the host I'm using). Run with
valgrind included below.

gdb --args ogrinfo /tmp/newdir/
...
(gdb) run
Starting program: /usr/local/bin/ogrinfo /tmp/newdir/
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
[New Thread 0x7fffe7757640 (LWP 988)]
[New Thread 0x7fffe6f56640 (LWP 989)]
[New Thread 0x7fffde755640 (LWP 990)]
[New Thread 0x7fffd5f54640 (LWP 991)]
[New Thread 0x7fffc5753640 (LWP 992)]
[New Thread 0x7fffc4f52640 (LWP 993)]
[New Thread 0x7fffb4751640 (LWP 994)]
[New Thread 0x7fffabf50640 (LWP 995)]
[New Thread 0x7fffab74f640 (LWP 996)]
[New Thread 0x7fffa2f4e640 (LWP 997)]
[New Thread 0x7fff9a74d640 (LWP 998)]
[New Thread 0x7fff91f4c640 (LWP 999)]
[New Thread 0x7fff8974b640 (LWP 1000)]
[New Thread 0x7fff78f4a640 (LWP 1001)]
[New Thread 0x7fff78749640 (LWP 1002)]
[New Thread 0x7fff6f5ff640 (LWP 1003)]

Thread 1 "ogrinfo" received signal SIGILL, Illegal instruction.
0x73773c9e in tiledb::common::ThreadPool::ThreadPool(unsigned long)
() from /lib/x86_64-linux-gnu/libtiledb.so.2.16




valgrind -s ogrinfo /tmp/newdir
==704== Memcheck, a memory error detector
==704== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
==704== Using Valgrind-3.18.1 and LibVEX; rerun with -h for copyright info
==704== Command: ogrinfo /tmp/newdir
==704==
INFO: Open of `/tmp/newdir'
  using driver `ESRI Shapefile' successful.
1: geo_melbourne (Polygon)
==704==
==704== HEAP SUMMARY:
==704== in use at exit: 25,486 bytes in 216 blocks
==704==   total heap usage: 15,761 allocs, 15,545 frees, 2,390,169 bytes
allocated
==704==
==704== LEAK SUMMARY:
==704==definitely lost: 0 bytes in 0 blocks
==704==indirectly lost: 0 bytes in 0 blocks
==704==  possibly lost: 544 bytes in 1 blocks
==704==still reachable: 24,942 bytes in 215 blocks
==704== suppressed: 0 bytes in 0 blocks
==704== Rerun with --leak-check=full to see details of leaked memory
==704==
==704== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)


ogrinfo /tmp/newdir
Illegal instruction (core dumped)

Cheers, Mike




On Sat, Feb 3, 2024 at 12:46 PM Even Rouault 
wrote:

> Michael,
>
> I'm wondering if there not might be something wrong with your build or
> runtime environment. Or there's something subtle, because that works fine
> for me with my dev build or in the ghcr.io/osgeo/gdal:alpine-normal-3.8.3
> Docker image
>
> Try running "valgrind ogrinfo /tmp/newdir/" or "gdb --args ogrinfo
> /tmp/newdir/" (type "run") to get more useful information
>
> Even
> Le 03/02/2024 à 02:35, Michael Sumner via gdal-dev a écrit :
>
> I'm getting Illegal instruction / core dumped on ogrinfo of a directory:
>
> ogr2ogr /tmp/newdir
> https://github.com/SymbolixAU/geojsonsf/raw/master/inst/examples/geo_melbourne.geojson
> -f "ESRI Shapefile"
>
> ogrinfo /tmp/newdir/
> Illegal instruction (core dumped)
>
> I've worked back through some docker images and it wasn't a problem in
> 3.6.0, but I'm getting it since 3.7.0 - or I'm doing something wrong
> entirely.
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] core dump on dir info

[Michael Sumner via gdal-dev]

I'm getting Illegal instruction / core dumped on ogrinfo of a directory:

ogr2ogr /tmp/newdir
https://github.com/SymbolixAU/geojsonsf/raw/master/inst/examples/geo_melbourne.geojson
-f "ESRI Shapefile"

ogrinfo /tmp/newdir/
Illegal instruction (core dumped)

I've worked back through some docker images and it wasn't a problem in
3.6.0, but I'm getting it since 3.7.0 - or I'm doing something wrong
entirely.

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Virtual Raster Tile Index (VRTTI) driver, and associated gdaltindex improvements

[Michael Sumner via gdal-dev]

just to follow up, I got it all working in latest GDAL:

script at

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

creates dsn

/vsicurl/https://github.com/mdsumner/cog-example/raw/main/gti/cop90.gti.fgb

which works nicely, thanks!

Cheers, Mike



On Tue, Jan 30, 2024 at 10:57 PM Even Rouault 
wrote:

>
> as you could also attach the XRES, YRES, BANDCOUNT, etc. metadata items
> that would help the GTI driver avoid probing one of the tiles.
>
> PS: OGRLayer::SetMetadataItem(key, value) only works on FlatGeoBuf and
> GPKG drivers with latest master as well
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Virtual Raster Tile Index (VRTTI) driver, and associated gdaltindex improvements

[Michael Sumner via gdal-dev]

ah thanks, all very helpful - no it's not done with master (for various
reasons), I'll follow up with details if relevant, mostly I was just
excited to get a useable workflow for the entire process.

Cheers, Mike

On Tue, 30 Jan 2024, 22:46 Even Rouault,  wrote:

> 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 "", line 1, in 
>   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 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

Re: [gdal-dev] Virtual Raster Tile Index (VRTTI) driver, and associated gdaltindex improvements

[Michael Sumner via gdal-dev]

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 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Layerstack multiple images, of different extents and variable number of bands

[Michael Sumner via gdal-dev]

I would start with

gdalwarp  out.tif

set -ts to something small at first to get a visual like -ts 1024 0   (the
y zero means the aspect ratio is figured out sanely from the x size  - with
no arguments you get best-preserving grid from all resolved inputs - set
some or 3 of -t_srs -tr -te -ts for a specific grid spec).

Just make sure the file with the largest number of bands is first in the
input list (separated by spaces), if actual data overlaps in matching
sequential bands last-in wins (IIUC).

Cheers, Mike









On Wed, Jan 24, 2024 at 4:02 AM lefsky--- via gdal-dev <
gdal-dev@lists.osgeo.org> wrote:

> I've been incorporating gdal into my processing pipelines and have a
> question about a feature.
>
> I need to be able to layer stack multiple images with different extents
> and varying numbers of bands, usually exceeding 1. I haven't been able to
> find a single command to do this in gdal. Am I missing something? If there
> is no such feature, has anyone developed a script out there that can
> robustly do all of those things (i.e. with checking for violation of
> assumptions about files and error handling)?
>
> Michael
>
> --
> Michael Lefsky (He/His)
> Home Location: HVHF+GH
> Cell: 970-980-9036
> http://www.researcherid.com/rid/A-7224-2009
>
> *“for being prematurely, and worse, intuitively right — there’s a heavy
> price. But for being wrong — no, not so long as you’re wrong in a pack."
> Gary Brecher / Portis*
>
> *I acknowledge that I live and work on stolen land. This is the land of
> the Cheyenne, Arapaho, Ute, and Ocheithi Sakowin people. To learn more
> about these nations, please visit;
> http://www.utemountainutetribe.com/
> http://www.cheyennenation.com/
> https://cheyenneandarapaho-nsn.gov/
> https://native-land.ca/
>
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] help please with close/reopen and GDALSubdatasetInfo

[Michael Sumner via gdal-dev]

Awesome, thank you!  I was a bit lost - learnt a lot from exploring this
though.

Cheers, Mike



On Mon, Dec 4, 2023 at 12:11 AM Even Rouault 
wrote:

> Hi Michael,
>
> It was most missing a "break;" statement when the subdataset of interest
> has been found to exit the loop. Without the break, the loop will
> continue to read papszSubdatasets[] which has been invalidated by the
> poSrcDS->ReleaseRef(). Ah, memory unsafe languages :-). You also had a
> few memory leaks.
>
> Fixed version:
>
> #include "gdal.h"
> #include "gdal_priv.h"
> #include 
> int main(int argc, char **argv) {
>GDALAllRegister();
>
>auto poSrcDS =
>  GDALDataset::Open(argv[1], GDAL_OF_RASTER, nullptr, nullptr, nullptr);
>if (poSrcDS == nullptr)
>{
>  return 0;
>}
>char **papszSubdatasets = poSrcDS->GetMetadata("SUBDATASETS");
>int nSubdatasets = CSLCount(papszSubdatasets);
>if (nSubdatasets > 0)
>{
>  for (int j = 0; j < nSubdatasets; j += 2)
>  {
>char* pszSubdatasetSource = CPLStrdup(strstr(papszSubdatasets[j],
> "=") + 1);
>GDALSubdatasetInfoH info =
> GDALGetSubdatasetInfo(pszSubdatasetSource);
>char* component = info ?
> GDALSubdatasetInfoGetSubdatasetComponent(info) : NULL;
>const bool bFound = component && EQUAL(argv[2], component);
>CPLFree(component);
>GDALDestroySubdatasetInfo(info);
>if ( bFound) {
>  std::cout << pszSubdatasetSource << "\n";
>  poSrcDS->ReleaseRef();
>  poSrcDS = GDALDataset::Open(pszSubdatasetSource,
> GDAL_OF_RASTER, nullptr, nullptr, nullptr);
>  CPLFree(pszSubdatasetSource);
>  break;
>}
>else {
>  CPLFree(pszSubdatasetSource);
>}
>  }
>}
>
>poSrcDS->ReleaseRef();
>return 1;
> }
>
> Even
>
>
> Le 03/12/2023 à 05:10, Michael Sumner via gdal-dev a écrit :
> > #include "gdal.h"
> > #include "gdal_priv.h"
> > #include 
> > int main(int argc, char **argv) {
> >   GDALAllRegister();
> >
> >   auto poSrcDS =
> > GDALDataset::Open(argv[1], GDAL_OF_RASTER, nullptr, nullptr,
> nullptr);
> >   if (poSrcDS == nullptr)
> >   {
> > return 0;
> >   }
> >   char **papszSubdatasets = poSrcDS->GetMetadata("SUBDATASETS");
> >   int nSubdatasets = CSLCount(papszSubdatasets);
> >   char *pszSubdatasetSource = nullptr;
> >   if (nSubdatasets > 0)
> >   {
> > for (int j = 0; j < nSubdatasets; j += 2)
> > {
> >   pszSubdatasetSource = CPLStrdup(strstr(papszSubdatasets[j], "=")
> > + 1);
> >   GDALSubdatasetInfoH info =
> > GDALGetSubdatasetInfo(pszSubdatasetSource);
> >   if ( EQUAL(argv[2],
> > GDALSubdatasetInfoGetSubdatasetComponent(info))) {
> > std::cout << pszSubdatasetSource << "\n";
> > poSrcDS->ReleaseRef();
> > poSrcDS = GDALDataset::Open(pszSubdatasetSource,
> > GDAL_OF_RASTER, nullptr, nullptr, nullptr);
> > CPLFree(pszSubdatasetSource);
> > GDALDestroySubdatasetInfo(info);
> >
> >   }
> > }
> >   }
> >
> >   poSrcDS->ReleaseRef();
> >   return 1;
> > }
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] help please with close/reopen and GDALSubdatasetInfo

[Michael Sumner via gdal-dev]

May I please ask for assistance with this code?  I'm trying to close a
dataset with subdatasets and then reopen if the requested subdataset is
found.

It segfaults for all but the last subdataset name ... and I'm at a loss for
what I'm doing wrong.

Thank you.

Code below, and this gist documents it with a little more details:
https://gist.github.com/mdsumner/2ffc9b4746454fa0152b4a250900e535

Cheers, Mike


#include "gdal.h"
#include "gdal_priv.h"
#include 
int main(int argc, char **argv) {
  GDALAllRegister();

  auto poSrcDS =
GDALDataset::Open(argv[1], GDAL_OF_RASTER, nullptr, nullptr, nullptr);
  if (poSrcDS == nullptr)
  {
return 0;
  }
  char **papszSubdatasets = poSrcDS->GetMetadata("SUBDATASETS");
  int nSubdatasets = CSLCount(papszSubdatasets);
  char *pszSubdatasetSource = nullptr;
  if (nSubdatasets > 0)
  {
for (int j = 0; j < nSubdatasets; j += 2)
{
  pszSubdatasetSource = CPLStrdup(strstr(papszSubdatasets[j], "=") + 1);
  GDALSubdatasetInfoH info = GDALGetSubdatasetInfo(pszSubdatasetSource);
  if ( EQUAL(argv[2], GDALSubdatasetInfoGetSubdatasetComponent(info))) {
std::cout << pszSubdatasetSource << "\n";
poSrcDS->ReleaseRef();
poSrcDS = GDALDataset::Open(pszSubdatasetSource, GDAL_OF_RASTER,
nullptr, nullptr, nullptr);
CPLFree(pszSubdatasetSource);
GDALDestroySubdatasetInfo(info);

      }
}
  }

  poSrcDS->ReleaseRef();
  return 1;
}

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] slow translate with OVERVIEW_LEVEL=NONE when no overviews exist

[Michael Sumner via gdal-dev]

When I translate this GeoTIFF to 10% original size, it's very very slow if
OVERVIEW_LEVEL=NONE is set.

The GeoTIFF has no overviews.

export dsn="/vsicurl/
https://github.com/mdsumner/cog-example/raw/main/cog/sentinel-image.tif;
## takes *forever*
gdal_translate $dsn  out.tif -outsize 219 226 -oo OVERVIEW_LEVEL=NONE

## works fast as expected
gdal_translate $dsn  out.tif -outsize 219 226

I found this in a package that has a global open option setting for
overviews for subsampling with RasterIO(). There we could detect that
overviews are not present, or even just not set it at all ... (better to
use GDALwarp() for general cases IMO).

But, is there a good reason why setting that open option affects translate
for sources with no overviews?

I'm in master at ebac6b74a429fa6eeeb94edd074b6cb60072a35f

Cheers, Mike



-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] format for topology?

[Michael Sumner via gdal-dev]

Are there any formats that record "coverage" topology? What I'm worried
about is when shapes are encoded as blob geometry with initially identical
 coordinates at shared vertices, is there any process that can validate or
record that particular coords should have the same values even after
reprojection or decomposition/rebuild?

I know how to do this in custom ways, I'm hoping to find some off-the-shelf
process to at least refer to (I'm a bit pressed for time so I'm asking
without doing the extra homework that's needed).  I would build an indexer
of the unique coordinates and how they relate to the blob geometries before
and after and do checks there, but I don't know of an existing process that
guarantees this.

I understand that precision-setting is used to prevent numeric loss, but is
there any formalization about that loss? (i.e.  different instances of a
vertex really do stay unique). I guess the result is going to be pretty
good if you choose the precision appropriately and modifications aren't
extreme?

Thanks for consideration!

Cheers, Mike






--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Using ogr2ogr with limited memory

[Michael Sumner via gdal-dev]

what's the task? what about batching the geometry and or fields? can you
run on just the first feature, does that work?

how many features, how big is the task?

On Wed, 27 Sept 2023, 13:16 Scott via gdal-dev, 
wrote:

> Any tips for using ogr2ogr to use only a specified amount of RAM? I'm
> not having any luck.
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] resampling algs warp vs translate

[Michael Sumner via gdal-dev]

Can we aspire to translate having the same set of resampling algorithms as
the warper?

I see the warper adds min, max, mod, q1, q3, sum

I especially wanted sum for OVERVIEW_RESAMPLING in COG, and I can see where
it's done and  ...  can maybe see my way through that ... but the 600 lines
of code in GDALResampleChunk_AverageOrRMS_T kind of reins in my confidence.

https://github.com/OSGeo/gdal/blob/26e7ff91504e0303434f440e6d19e58f2a206b45/gcore/overview.cpp#L1211C15-L1211C47

It came up as a need for statistical properties required for
overview-publishers, and while I think we can incrementally build up the
SUM-pyramid with external overivews via CLI it would be awesome to have the
same stats for the two resampling regimes.

Thank you



-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] adding metadata -mo in a domain?

[Michael Sumner]

cool thanks, I'll deal with that 



On Thu, 21 Sept 2023, 18:48 Even Rouault, 
wrote:

>
> >
> > which don't seem to exist. Is that just an outdated comment?
>
> yes, as far as I can see from history those functions have never
> existed. This comment dates back to the initial commit and they didn't
> even exist at that time.
>
>
> --
> 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

Re: [gdal-dev] adding metadata -mo in a domain?

[Michael Sumner]

I see this note about

GDALTranslateOptionsSetMetadataOptions()

GDALTranslateOptionsAddMetadataOptions()

which don't seem to exist. Is that just an outdated comment?

https://github.com/OSGeo/gdal/blob/8122e2c39f72317a803cd24ee1fb5d6a44a97596/apps/gdal_translate_lib.cpp#L200-L201

Thanks!



On Wed, Sep 20, 2023 at 6:40 PM Even Rouault 
wrote:

> Sounds OK to me
> Le 20/09/2023 à 02:32, Michael Sumner a écrit :
>
> Thanks Even,  could add a  AttachDomainMetadata here that expected "-dmo
> DOMAIN:META-TAG=VALUE"
>
>
> https://github.com/OSGeo/gdal/blob/a4b57c22606be2258196a933d310f301e3ece93f/apps/gdal_translate_lib.cpp#L2625
>
> I'll try that, or is it too much over-reach?
>
> Ultimately it's for vrt:// con so could do
>
> vrt://{dsn}?dmo=DOMAIN:META-TAG=VALUE
>
> Cheers, Mike
>
>
>
> On Wed, Sep 20, 2023 at 10:19 AM Even Rouault 
> wrote:
>
>> Not from cli. From the API, it is dataset.SetMetadataItem(key, value,
>> domain)
>> Le 20/09/2023 à 02:14, Michael Sumner a écrit :
>>
>> When we add metadata with -mo can we somehow specify the DOMAIN it would
>> be added to, or be instantiated with?
>>
>> gdal_translate gcore/data/float32.tif /tmp/add_md.vrt -mo PHONY=important
>>
>> I would like to tweak that so that the item was in a new domain, i.e.
>>
>>
>>   
>> Area
>>   
>>   
>> BAND
>>   
>>   
>>  important
>>   
>>
>> Is that possible from cli?
>>
>> Thank you
>>
>>
>> --
>> Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> Hobart, Australia
>> e-mail: mdsum...@gmail.com
>>
>> ___
>> gdal-dev mailing 
>> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> -- http://www.spatialys.com
>> My software is free, but my time generally not.
>>
>>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] adding metadata -mo in a domain?

[Michael Sumner]

Thanks Even,  could add a  AttachDomainMetadata here that expected "-dmo
DOMAIN:META-TAG=VALUE"

https://github.com/OSGeo/gdal/blob/a4b57c22606be2258196a933d310f301e3ece93f/apps/gdal_translate_lib.cpp#L2625

I'll try that, or is it too much over-reach?

Ultimately it's for vrt:// con so could do

vrt://{dsn}?dmo=DOMAIN:META-TAG=VALUE

Cheers, Mike



On Wed, Sep 20, 2023 at 10:19 AM Even Rouault 
wrote:

> Not from cli. From the API, it is dataset.SetMetadataItem(key, value,
> domain)
> Le 20/09/2023 à 02:14, Michael Sumner a écrit :
>
> When we add metadata with -mo can we somehow specify the DOMAIN it would
> be added to, or be instantiated with?
>
> gdal_translate gcore/data/float32.tif /tmp/add_md.vrt -mo PHONY=important
>
> I would like to tweak that so that the item was in a new domain, i.e.
>
>
>   
> Area
>   
>   
> BAND
>   
>   
>  important
>   
>
> Is that possible from cli?
>
> Thank you
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] adding metadata -mo in a domain?

[Michael Sumner]

When we add metadata with -mo can we somehow specify the DOMAIN it would be
added to, or be instantiated with?

gdal_translate gcore/data/float32.tif /tmp/add_md.vrt -mo PHONY=important

I would like to tweak that so that the item was in a new domain, i.e.


  
Area
  
  
BAND
  
  
 important
  

Is that possible from cli?

Thank you


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] -scale default

[Michael Sumner]

nice, thanks so much for the extended description.

Cheers, Mike

On Fri, 8 Sept 2023, 21:34 Even Rouault,  wrote:

> Michael,
>
> It turns out that gdal_translate was using 255.999 since the origins in
> 2002 :
> https://github.com/OSGeo/gdal/commit/47d859efb25da249b5d62eb8619c36134ae76572#diff-d57b4553312c88a736821448adc155e6797653b0d902bd42f9310572aa71f0b0R215
>
> I suspect the reasoning was that, if you rounded values down (floor) and
> used 255 as the scaleDstMax, and you remapped [0,65535] to [0,255] for
> example, then 65534. / (65535. / 255) = 254.996.. would be rounded to 254.
> So you would only get 255 as the output for 65535 as the input, which seems
> "unfair". But GDAL actually rounds to closest, so if using dstMax = 65535,
> all input values in range [65535 - 128, 65535] get scaled to 255, which is
> OK to me. Of course that last "bucket" and the first one at 0 (input values
> in [0,128]) are half-width of the ones in the middle (input values in
> [129,129+256] get scaled to 1) but I don't think we can do much about that
> (discrete values vs intervals)
>
> One downside of using dstMax=255.999 is that if your input range is [0,
> 255], then -scale was not the identity (254 got transformed to 253), which
> is kind of an issue! Hence I've changed dstMax to be 255 in
> https://github.com/OSGeo/gdal/pull/8367
>
> I've also noticed during testing a more subtle consistency issue when if
> you did for example "gdal_translate autotest/gcore/data/byte.tif out.asc
> -scale 0 1 0 1.5 -of aaigrid", you would get values outside of the [0, 255]
> range because the VRT machinery was ignoring the constraints of the
> VRTRasterBand data type (Byte), and was only honouring the buffer data type
> of the RasterIO() request from the AAIGrid driver (Int32, always used when
> translating from source bands having integer data types, but with the
> reasoning that if you query a Byte band then you would get always value in
> [0,255] range)
>
> Even
> Le 08/09/2023 à 04:20, Michael Sumner a écrit :
>
> Hi, I'm expecting `gdal_translate -scale` to emit values in the range
> 0,255 but it seems to be targeting 0,256.  (All works as expected when
> using explicit src_min src_max dst_min dst_max).
>
> To see the output range on a simple example (see code at link below in
> case email garbles):
>
> from osgeo import gdal
> ds = gdal.Translate("/vsimem/scl.tif", "autotest/gcore/data/float32.tif",
> options = "-scale")
> ds.GetRasterBand(1).ComputeRasterMinMax()
> ## (0.0, 255.99899291992188)
>
> Is that expected, am I missing something? I've tried variations on the
> input values and output types.
>
> https://gist.github.com/mdsumner/ee4103d8616b9aa341e82c46b44a8c8c
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Research Software Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 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
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] -scale default

[Michael Sumner]

Hi, I'm expecting `gdal_translate -scale` to emit values in the range 0,255
but it seems to be targeting 0,256.  (All works as expected when using
explicit src_min src_max dst_min dst_max).

To see the output range on a simple example (see code at link below in case
email garbles):

from osgeo import gdal
ds = gdal.Translate("/vsimem/scl.tif", "autotest/gcore/data/float32.tif",
options = "-scale")
ds.GetRasterBand(1).ComputeRasterMinMax()
## (0.0, 255.99899291992188)

Is that expected, am I missing something? I've tried variations on the
input values and output types.

https://gist.github.com/mdsumner/ee4103d8616b9aa341e82c46b44a8c8c

Cheers, Mike


-- 
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] multidim arraySpecs/GetView

[Michael Sumner]

Excellent, thanks so much, confirmed it's working.

Cheers, Mike

On Fri, Sep 1, 2023 at 6:14 AM Even Rouault 
wrote:

> Michael,
>
> Hi, hoping to check my syntax ... or expectations ... I'm expecting this
> to subset by step the 20x20 array in this netcdf, but getting no change in
> python MultiDimTranslate or at command line:
>
> ds = gdal.OpenEx('gdrivers/data/netcdf/byte_no_cf.nc',
> gdal.OF_MULTIDIM_RASTER)
>
> nds = gdal.MultiDimTranslate("/vsimem/array_view.zarr", ds, format="Zarr",
> arraySpecs=['name=Band1,view=[::2,::4]'])
> nds.GetRootGroup().OpenMDArray("Band1").GetShape()
> # (20, 20)
>
> This was a bug. Fixed per https://github.com/OSGeo/gdal/pull/8297
>
> A way of working it around is to use  scaleAxesSpecs=['x(4)', 'y(2)']
> instead
>
> Even
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] multidim arraySpecs/GetView

[Michael Sumner]

Hi, hoping to check my syntax ... or expectations ... I'm expecting this to
subset by step the 20x20 array in this netcdf, but getting no change in
python MultiDimTranslate or at command line:

ds = gdal.OpenEx('gdrivers/data/netcdf/byte_no_cf.nc',
gdal.OF_MULTIDIM_RASTER)

nds = gdal.MultiDimTranslate("/vsimem/array_view.zarr", ds, format="Zarr",
arraySpecs=['name=Band1,view=[::2,::4]'])
nds.GetRootGroup().OpenMDArray("Band1").GetShape()
# (20, 20)

Directly with the python GetView api it works

ds.GetRootGroup().OpenMDArray("Band1").GetView("[::2,::4]").GetShape()
# (10, 5)

Am I specifying the view wrong in MultiDimTranslate?

Thanks, Mike

I've put fuller notes here in case the email is a problem:
https://gist.github.com/mdsumner/028d089ad9a960267347f891c058fc9a

--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] convert NETCDF to Geotiff upside down

[Michael Sumner]

just a follow up, as far as I can tell the lon,lat arrays are actually
natively in global Mercator, so a four-number extent would suffice  but
these arrays have been materialized rather than preserve that compact
origin (it's very common sadly). You might feed that back to the providers
and/or other users of the data, shipping these lon,lat arrays is
unfortunate when simple georeferencing (extent, expressed as affine
transform and a crs) would do. There's more metadata than data in these
files ...

I've used many tricks to workaround bad formats like these, but now I'd
just say to shove it through the warper because that is now so fast and
reliable.

Cheers, MIke



On Thu, Aug 10, 2023 at 7:17 AM Michael Sumner  wrote:

> there is no georeferencing with this one, but as you say it has internal
> coordinates (which GDAL interprets as "geolocation arrays")
>
> so, run it through the warper i.e.
>
> gdalwarp NETCDF:GLOBCOMPLIR_nc.2023080819:data ofn.tif
>
> With that you can (and should) set your desired extent (-te), dimensions
> (-ts) [or resolution -tr], and crs (-t_crs) but GDAL will produce an output
> that works by heuristics if you specify none (or some) of those.
>
> I did the following just for fun to see if the non-downloading pathway
> would work it (not sure if the config is required or not, for testing
> purposes I set dimension to 1024,0 which means GDAL figures out a sensible
> aspect ratio - just remove -ts when you're happy with the result for full
> resolution. )
>
> gdalwarp --config AWS_NO_SIGN_REQUEST "YES"
> NETCDF:"/vsis3/noaa-gmgsi-pds/GMGSI_LW/2023/08/08/19/GLOBCOMPLIR_nc.2023080819":data
> -ts 1024 0 GLOBCOMPLIR_nc.2023080819_test.tif
>
> (some details might not work depending on your GDAL version)
>
> Cheers, Mike
>
> On Thu, Aug 10, 2023 at 6:30 AM Dori  wrote:
>
>> Hi,
>> I tried to convert the NETCDF stored in
>> https://noaa-gmgsi-pds.s3.amazonaws.com/index.html#GMGSI_LW/202308/08/19/GLOBCOMPLIR_nc.2023080819
>>  to Geotiff format using following command:
>>
>> gdal_translate 0ol GTiff -a_srs EPSG:4326 -aullr -180 -72.7 179.9 -72.7
>> -co “COMPRESS=LZW” NETCDF:GLOBCOMPLIR_nc.2023080819:data ofn.tif
>>
>> ofn.tif is upside down.  North latitude points are in South and viceversa.
>>
>> How can I solve this upside down output?
>> What causes this?
>>
>> I did notice that the points in the NETCDF file are not equidistant  what
>> could be my issue.
>>
>> Thank you for any hint on how to go about getting the correct Geotiff
>> output
>>
>> ___
>> 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
>


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] convert NETCDF to Geotiff upside down

[Michael Sumner]

there is no georeferencing with this one, but as you say it has internal
coordinates (which GDAL interprets as "geolocation arrays")

so, run it through the warper i.e.

gdalwarp NETCDF:GLOBCOMPLIR_nc.2023080819:data ofn.tif

With that you can (and should) set your desired extent (-te), dimensions
(-ts) [or resolution -tr], and crs (-t_crs) but GDAL will produce an output
that works by heuristics if you specify none (or some) of those.

I did the following just for fun to see if the non-downloading pathway
would work it (not sure if the config is required or not, for testing
purposes I set dimension to 1024,0 which means GDAL figures out a sensible
aspect ratio - just remove -ts when you're happy with the result for full
resolution. )

gdalwarp --config AWS_NO_SIGN_REQUEST "YES"
NETCDF:"/vsis3/noaa-gmgsi-pds/GMGSI_LW/2023/08/08/19/GLOBCOMPLIR_nc.2023080819":data
-ts 1024 0 GLOBCOMPLIR_nc.2023080819_test.tif

(some details might not work depending on your GDAL version)

Cheers, Mike

On Thu, Aug 10, 2023 at 6:30 AM Dori  wrote:

> Hi,
> I tried to convert the NETCDF stored in
> https://noaa-gmgsi-pds.s3.amazonaws.com/index.html#GMGSI_LW/202308/08/19/GLOBCOMPLIR_nc.2023080819
>  to Geotiff format using following command:
>
> gdal_translate 0ol GTiff -a_srs EPSG:4326 -aullr -180 -72.7 179.9 -72.7
> -co “COMPRESS=LZW” NETCDF:GLOBCOMPLIR_nc.2023080819:data ofn.tif
>
> ofn.tif is upside down.  North latitude points are in South and viceversa.
>
> How can I solve this upside down output?
> What causes this?
>
> I did notice that the points in the NETCDF file are not equidistant  what
> could be my issue.
>
> Thank you for any hint on how to go about getting the correct Geotiff
> output
>
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] STACIT fail to open s3://..tif

[Michael Sumner]

Ok so I just did this in staticdataset.cpp and it works ... I'll follow up
on github.

else if (STARTS_WITH(osFilename.c_str(), "s3://"))
{
osRet = "/vsis3/";
osRet += osFilename.substr(strlen("s3://"));
}


Thanks!



On Wed, Aug 2, 2023 at 5:23 PM Michael Sumner  wrote:

> I'm using a derived subdataset dsn to try to open a STAC query:
>
> gdalinfo "STACIT:\"
> https://explorer.sandbox.dea.ga.gov.au/stac/search?bbox=122.1,-18.28,122.48,-17.91=2020-01-01/2020-02-29\
> ":asset=level3,collections=ga_ls_landcover_class_cyear_2"
>
> ERROR 1: Cannot open
> s3://dea-public-data/derivative/ga_ls_landcover_class_cyear_2/1-0-0/2020/x_-11/y_-20/ga_ls_landcover_class_cyear_2_1-0-0_au_x-11y-20_2020-01-01_level3.tif
> to retrieve band characteristics
> gdalinfo failed - unable to open 'STACIT:"
> https://explorer.sandbox.dea.ga.gov.au/stac/search?bbox=122.1,-18.28,122.48,-17.91=2020-01-01/2020-02-29
> ":asset=level3,collections=ga_ls_landcover_class_cyear_2'.
>
> Am I doing something wrong with the STAC approach?  Similar tasks are
> working well with Microsoft Planetary Computer (but I'm generally out of my
> depth with STAC itself).
>
> The .tif itself seems fine via the /vsis3/ syntax
>
> gdalinfo
> /vsis3/dea-public-data/derivative/ga_ls_landcover_class_cyear_2/1-0-0/2020/x_-11/y_-20/ga_ls_landcover_class_cyear_2_1-0-0_au_x-11y-20_2020-01-01_level3.tif
>
>
> I've put the calls here in case email problems with the strings:
> https://gist.github.com/mdsumner/90f5b509c5cb32bc0f1a26c695baa0bf
>
> I'm using gdal dev GDAL 3.8.0dev-bba8edfc13, released 2023/08/01 (debug
> build)
>
> Thanks, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] STACIT fail to open s3://..tif

[Michael Sumner]

I'm using a derived subdataset dsn to try to open a STAC query:

gdalinfo "STACIT:\"
https://explorer.sandbox.dea.ga.gov.au/stac/search?bbox=122.1,-18.28,122.48,-17.91=2020-01-01/2020-02-29\
":asset=level3,collections=ga_ls_landcover_class_cyear_2"

ERROR 1: Cannot open
s3://dea-public-data/derivative/ga_ls_landcover_class_cyear_2/1-0-0/2020/x_-11/y_-20/ga_ls_landcover_class_cyear_2_1-0-0_au_x-11y-20_2020-01-01_level3.tif
to retrieve band characteristics
gdalinfo failed - unable to open 'STACIT:"
https://explorer.sandbox.dea.ga.gov.au/stac/search?bbox=122.1,-18.28,122.48,-17.91=2020-01-01/2020-02-29
":asset=level3,collections=ga_ls_landcover_class_cyear_2'.

Am I doing something wrong with the STAC approach?  Similar tasks are
working well with Microsoft Planetary Computer (but I'm generally out of my
depth with STAC itself).

The .tif itself seems fine via the /vsis3/ syntax

gdalinfo
/vsis3/dea-public-data/derivative/ga_ls_landcover_class_cyear_2/1-0-0/2020/x_-11/y_-20/ga_ls_landcover_class_cyear_2_1-0-0_au_x-11y-20_2020-01-01_level3.tif


I've put the calls here in case email problems with the strings:
https://gist.github.com/mdsumner/90f5b509c5cb32bc0f1a26c695baa0bf

I'm using gdal dev GDAL 3.8.0dev-bba8edfc13, released 2023/08/01 (debug
build)

Thanks, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] configuration option GDAL_RASTERIO_RESAMPLING

[Michael Sumner]

Ah great thanks, interested to help flush these out.

Cheers, Mike






On Sun, 28 May 2023, 10:34 Daniel Baston,  wrote:

> Hi Mike,
>
> The content on the page you linked is for the most part a simple transfer
> of what was available on the wiki, plus a few options that were previously
> defined on driver or other pages. It is not complete.
>
> There are also about 500 configuration options that are used in the code
> via CPLGetConfigOption that do not appear in the documentation,
> GDAL_RASTERIO_RESAMPLING being one of them. Some of these are clearly not
> intended to be part of the GDAL interface, but there's still quite a bit
> that can be added to the documentation. I can create an issue and post the
> script I used to scan for these if others are interested in filling some of
> the gaps.
>
> Dan
>
>
>
>
>
>
>
> On Sat, May 27, 2023 at 8:18 PM Michael Sumner  wrote:
>
>> Hi, I'm very happy to see the refactor of documentation for configuration
>> and environment options, primarily by Dan Baston in recent weeks - thanks!
>>
>> https://gdal.org/user/configoptions.html
>>
>> Does the GDAL_RASTERIO_RESAMPLING option belong in this list?  I'm
>> exploring its use versus setting resampling with the API, and I gather also
>> that there is not an analogous config option for the warper.
>>
>> Thank you, Mike
>>
>> --
>> Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> Hobart, Australia
>> e-mail: mdsum...@gmail.com
>> ___
>> gdal-dev mailing list
>> gdal-dev@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] configuration option GDAL_RASTERIO_RESAMPLING

[Michael Sumner]

Hi, I'm very happy to see the refactor of documentation for configuration
and environment options, primarily by Dan Baston in recent weeks - thanks!

https://gdal.org/user/configoptions.html

Does the GDAL_RASTERIO_RESAMPLING option belong in this list?  I'm
exploring its use versus setting resampling with the API, and I gather also
that there is not an analogous config option for the warper.

Thank you, Mike

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] request limit on netcdf server ?

[Michael Sumner]

I'm having trouble with remote dsn from here:

https://esgf.nci.org.au/thredds/catalog/esgcet/cmip5/output1/CSIRO-BOM/ACCESS1-3/rcp85/day/ocean/day/r1i1p1/cmip5.output1.CSIRO-BOM.ACCESS1-3.rcp85.day.ocean.day.r1i1p1.v20120413.html#cmip5.output1.CSIRO-BOM.ACCESS1-3.rcp85.day.ocean.day.r1i1p1.v20120413

The support tells me I have been blocked for requests from my IP that
exceed 300 in 5min. Apparently, daily you are cleared from the blocklist
but with experiments today I've had two machines blocked (and now manually
cleared). Still, it's hard to test :0

I've stripped down to this dsn, I thought limiting  the driver might limit
the requests because it won't cascade down the driver identification (but
I'm not sure if -if adds anymore than the "NETCDF:" is already doing). At
any rate I need to target a subdataset so I need that syntax.

"vrt://NETCDF:\"/vsicurl/
https://esgf.nci.org.au/thredds/fileServer/master/CMIP5/output1/CSIRO-BOM/ACCESS1-3/rcp85/day/ocean/day/r1i1p1/v20120413/tos/tos_day_ACCESS1-3_rcp85_r1i1p1_20560101-20651231.nc\
":tos?bands=1=NETCDF"

It takes quite a while to read one band, but the warper and RasterIO result
is good.  I think the machine is now blocked again, I wasn't able to also
run gdalinfo after warper and straight read.

I'm just doing standalone from-scratch opens, not with a persistent dataset
- I assume that would reduce the request load, which I'll try in python
when I have access again.

Questions:

1. Are there other things I can pursue to reduce my request load on the
server?

2. Can this actually be faster?  RasterIO for one band 360x300 was
minutes.  (I appreciate this might be better suited to mdim, but that seems
similarly problematic in terms of requests and performance).

I've tried the config options   GDAL_NETCDF_VERIFY_DIMS / STRICT and
GDAL_NETCDF_IGNORE_XY_AXIS_NAME_CHECK=/YES and that does seem to have
helped the warper request but RasterIO was still slow.
.
This was on GDAL 3.8.0dev-414270a94b-dirty, released 2023/05/09 (debug
build)

I see this repeated warning

Warning 1: dimension #2 (i) is not a Longitude/X dimension.
Warning 1: dimension #1 (j) is not a Latitude/Y dimension.
Warning 1: dimension #1 (i) is not a Longitude/X dimension.
Warning 1: dimension #0 (j) is not a Latitude/Y dimension.
Warning 1: dimension #1 (i) is not a Longitude/X dimension.
Warning 1: dimension #0 (j) is not a Latitude/Y dimension.

Thank you

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDALOpenShared: warning about "filename" vs "description"

[Michael Sumner]

Thank you.

Cheers, Mike



On Sun, Apr 30, 2023 at 1:07 AM Even Rouault 
wrote:

> Hi Michael,
>
> Fixed per https://github.com/OSGeo/gdal/pull/7659
>
> Even
> Le 24/04/2023 à 05:09, Michael Sumner a écrit :
>
>
> When invoking  info on this vrt connection
>
> gdalinfo "vrt://
> https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2?a_ullr=0,9661,7180,0;
>
>
> I see this warning
>
> Warning 6: A dataset opened by GDALOpenShared should have the same
> filename (
> https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2)
> and description (/vsimem/http_1/NSW041500.jp2)
>
> The warning seems unreasonable given that we are quite far from a
> "filename" for a dsn? There's no warning on the bare url or with /vsicurl:
>
> gdalinfo
> https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2
>
> gdalinfo  /vsicurl/
> https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2
>
> ## GDAL version and provenance
>
> Ubuntu 20.04
>
> GDAL 3.7.0dev-717dcc0eed, released 2023/04/22 (debug build)
>
>
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDALOpenShared: warning about "filename" vs "description"

[Michael Sumner]

When invoking  info on this vrt connection

gdalinfo "vrt://
https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2?a_ullr=0,9661,7180,0;


I see this warning

Warning 6: A dataset opened by GDALOpenShared should have the same filename
(
https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2)
and description (/vsimem/http_1/NSW041500.jp2)

The warning seems unreasonable given that we are quite far from a
"filename" for a dsn? There's no warning on the bare url or with /vsicurl:

gdalinfo
https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2

gdalinfo  /vsicurl/
https://herbariumnsw-pds.s3-ap-southeast-2.amazonaws.com/images/NSW041500.jp2

## GDAL version and provenance

Ubuntu 20.04

GDAL 3.7.0dev-717dcc0eed, released 2023/04/22 (debug build)




-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] recent change?, issue with OGRGeometryFactory in R

[Michael Sumner]

I figured it out, I'd moved to from-source builds to /usr on a machine that
had apt installs to /usr/local so the stale headers were in the way.

A good lesson for me (don't do what I'm doing it's fraught haha) and just a
note to clear the air.

Cheers, Mike



On Sat, 1 Apr 2023, 21:09 Michael Sumner,  wrote:

> ah, thanks so much - maybe I'm seeing partial build problems on one
> machine, but at any rate this gives me enough to explore and figure it out.
>
> Best, Mike
>
> On Sat, 1 Apr 2023, 21:03 Even Rouault, 
> wrote:
>
>> Michael,
>>
>> There have indeed been recent changes in const correctness of a few
>> methods returning or accepting a OGRSpatialReference* (cf
>> https://github.com/OSGeo/gdal/pull/7381), including
>> OGRGeometryFactory::createmFromWkt().
>>
>> This requires a rebuild of code depending on the GDAL C++ API, as the
>> GDAL ABI generally changes between feature versions.
>>
>> Even
>> Le 01/04/2023 à 07:24, Michael Sumner a écrit :
>>
>> Hi,
>>
>> I'm here to ask for help finding where this issue is. If I use an earlier
>> build it still works (these are R packages sf and terra that are well
>> tested against previous releases).
>>
>> Has something changed with the headers structure for OGRGeometryFactory
>> perhaps? (I'm out of my depth for how to track this down).   If it's an
>> issue on the R side I'll fix it there, but mostly I'm here to ask what has
>> changed on the GDAL side to result in this symbol problem (or whatever it
>> is I'm doing wrong).
>>
>> using latest build  24afc75d3c3f2c85621e8b8d83f54a73b1d3ad33
>>
>> I see
>>
>> > library(sf)
>> Error: package or namespace load failed for ‘sf’ in dyn.load(file,
>> DLLpath = DLLpath, ...):
>>  unable to load shared object
>> '/usr/local/lib/R/site-library/sf/libs/sf.so':
>>   /usr/local/lib/R/site-library/sf/libs/sf.so: undefined symbol:
>> _ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
>> > library(terra)
>> Error: package or namespace load failed for ‘terra’ in dyn.load(file,
>> DLLpath = DLLpath, ...):
>>  unable to load shared object
>> '/usr/local/lib/R/site-library/terra/libs/terra.so':
>>   /usr/local/lib/R/site-library/terra/libs/terra.so: undefined symbol:
>> _ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
>>
>> Using a prior build eca6f7bf17fe500ef4c0ba9793e93ea305bab1fc it works.
>>
>> The offending code in terra is this (I think)
>>
>>
>> https://github.com/rspatial/terra/blob/ddcf84ce901a269ac13118d650efac16f89a8602/src/read_ogr.cpp#L694
>>
>> In sf it is
>>
>>
>> https://github.com/r-spatial/sf/blob/main/src/gdal_read.cpp#L582
>> Thank you.
>>
>> Cheers, Mike
>>
>>
>> --
>> Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> Hobart, Australia
>> e-mail: mdsum...@gmail.com
>>
>> ___
>> gdal-dev mailing 
>> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> -- 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

Re: [gdal-dev] recent change?, issue with OGRGeometryFactory in R

[Michael Sumner]

ah, thanks so much - maybe I'm seeing partial build problems on one
machine, but at any rate this gives me enough to explore and figure it out.

Best, Mike

On Sat, 1 Apr 2023, 21:03 Even Rouault,  wrote:

> Michael,
>
> There have indeed been recent changes in const correctness of a few
> methods returning or accepting a OGRSpatialReference* (cf
> https://github.com/OSGeo/gdal/pull/7381), including
> OGRGeometryFactory::createmFromWkt().
>
> This requires a rebuild of code depending on the GDAL C++ API, as the GDAL
> ABI generally changes between feature versions.
>
> Even
> Le 01/04/2023 à 07:24, Michael Sumner a écrit :
>
> Hi,
>
> I'm here to ask for help finding where this issue is. If I use an earlier
> build it still works (these are R packages sf and terra that are well
> tested against previous releases).
>
> Has something changed with the headers structure for OGRGeometryFactory
> perhaps? (I'm out of my depth for how to track this down).   If it's an
> issue on the R side I'll fix it there, but mostly I'm here to ask what has
> changed on the GDAL side to result in this symbol problem (or whatever it
> is I'm doing wrong).
>
> using latest build  24afc75d3c3f2c85621e8b8d83f54a73b1d3ad33
>
> I see
>
> > library(sf)
> Error: package or namespace load failed for ‘sf’ in dyn.load(file, DLLpath
> = DLLpath, ...):
>  unable to load shared object
> '/usr/local/lib/R/site-library/sf/libs/sf.so':
>   /usr/local/lib/R/site-library/sf/libs/sf.so: undefined symbol:
> _ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
> > library(terra)
> Error: package or namespace load failed for ‘terra’ in dyn.load(file,
> DLLpath = DLLpath, ...):
>  unable to load shared object
> '/usr/local/lib/R/site-library/terra/libs/terra.so':
>   /usr/local/lib/R/site-library/terra/libs/terra.so: undefined symbol:
> _ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
>
> Using a prior build eca6f7bf17fe500ef4c0ba9793e93ea305bab1fc it works.
>
> The offending code in terra is this (I think)
>
>
> https://github.com/rspatial/terra/blob/ddcf84ce901a269ac13118d650efac16f89a8602/src/read_ogr.cpp#L694
>
> In sf it is
>
>
> https://github.com/r-spatial/sf/blob/main/src/gdal_read.cpp#L582
> Thank you.
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 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

Re: [gdal-dev] recent change?, issue with OGRGeometryFactory in R

[Michael Sumner]

Actually... this seems to be a problem on just one machine, I would have
thought they are the same but fwiw this seems not to be a GDAL or R package
problem.

I'll check again carefully so probably ignore this, thanks!

Cheers, Mike

On Sat, 1 Apr 2023, 16:24 Michael Sumner,  wrote:

> Hi,
>
> I'm here to ask for help finding where this issue is. If I use an earlier
> build it still works (these are R packages sf and terra that are well
> tested against previous releases).
>
> Has something changed with the headers structure for OGRGeometryFactory
> perhaps? (I'm out of my depth for how to track this down).   If it's an
> issue on the R side I'll fix it there, but mostly I'm here to ask what has
> changed on the GDAL side to result in this symbol problem (or whatever it
> is I'm doing wrong).
>
> using latest build  24afc75d3c3f2c85621e8b8d83f54a73b1d3ad33
>
> I see
>
> > library(sf)
> Error: package or namespace load failed for ‘sf’ in dyn.load(file, DLLpath
> = DLLpath, ...):
>  unable to load shared object
> '/usr/local/lib/R/site-library/sf/libs/sf.so':
>   /usr/local/lib/R/site-library/sf/libs/sf.so: undefined symbol:
> _ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
> > library(terra)
> Error: package or namespace load failed for ‘terra’ in dyn.load(file,
> DLLpath = DLLpath, ...):
>  unable to load shared object
> '/usr/local/lib/R/site-library/terra/libs/terra.so':
>   /usr/local/lib/R/site-library/terra/libs/terra.so: undefined symbol:
> _ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
>
> Using a prior build eca6f7bf17fe500ef4c0ba9793e93ea305bab1fc it works.
>
> The offending code in terra is this (I think)
>
>
> https://github.com/rspatial/terra/blob/ddcf84ce901a269ac13118d650efac16f89a8602/src/read_ogr.cpp#L694
>
> In sf it is
>
>
> https://github.com/r-spatial/sf/blob/main/src/gdal_read.cpp#L582
> Thank you.
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] recent change?, issue with OGRGeometryFactory in R

[Michael Sumner]

Hi,

I'm here to ask for help finding where this issue is. If I use an earlier
build it still works (these are R packages sf and terra that are well
tested against previous releases).

Has something changed with the headers structure for OGRGeometryFactory
perhaps? (I'm out of my depth for how to track this down).   If it's an
issue on the R side I'll fix it there, but mostly I'm here to ask what has
changed on the GDAL side to result in this symbol problem (or whatever it
is I'm doing wrong).

using latest build  24afc75d3c3f2c85621e8b8d83f54a73b1d3ad33

I see

> library(sf)
Error: package or namespace load failed for ‘sf’ in dyn.load(file, DLLpath
= DLLpath, ...):
 unable to load shared object '/usr/local/lib/R/site-library/sf/libs/sf.so':
  /usr/local/lib/R/site-library/sf/libs/sf.so: undefined symbol:
_ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry
> library(terra)
Error: package or namespace load failed for ‘terra’ in dyn.load(file,
DLLpath = DLLpath, ...):
 unable to load shared object
'/usr/local/lib/R/site-library/terra/libs/terra.so':
  /usr/local/lib/R/site-library/terra/libs/terra.so: undefined symbol:
_ZN18OGRGeometryFactory13createFromWktEPKcP19OGRSpatialReferencePP11OGRGeometry

Using a prior build eca6f7bf17fe500ef4c0ba9793e93ea305bab1fc it works.

The offending code in terra is this (I think)

https://github.com/rspatial/terra/blob/ddcf84ce901a269ac13118d650efac16f89a8602/src/read_ogr.cpp#L694

In sf it is


https://github.com/r-spatial/sf/blob/main/src/gdal_read.cpp#L582
Thank you.

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Off topic COG question

[Michael Sumner]

I saw this : https://github.com/mapbox/COGDumper

On Sat, 1 Apr 2023, 02:54 Travis Kirstine,  wrote:

> Sorry to post this unrelated question to the GDAL group but I was hoping
> someone could provide some recommendations.  Does anybody know of a way to
> efficiently publish a WTMS / TMS from a COG source.  I know I can use
> something like MapServer / Mapcache or GeoServer can use the cog as a
> raster source but this seems like additional overhead as the COG internal
> tiling and overview structure is already optimized.  I know of titiler but
> that is about it.  https://github.com/developmentseed/titiler
>
> Thanks
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] question about python installation

[Michael Sumner]

Ah, awesome - thanks as ever for such great detail.

I've been trying with PYTHONPATH and sys.path, and I think I was doing it
wrong ...   but the simplicity of the /usr prefix had not occurred to me,
I'll be doing that.

If I come back around with further insights I'll share them, but this is
working for me now.

cmake .. -DCMAKE_INSTALL_PREFIX=/usr

Cheers, Mike





On Wed, Mar 8, 2023 at 1:35 PM Even Rouault 
wrote:

> Michael,
>
> My experience is that understanding how and where to install Python stuff
> is utterly difficult. And this is even more complicated on Debian/Ubuntu
> that have many patches regarding sys.path and distutils, plus the fact that
> distutils is deprecated and going to be removed in Python 3.12, without a
> clear replacement. Just look at all the horrible logic in
> swig/python/CMakeLists.txt from line 428 to 504 +
> swig/python/trimmedsysconfig.py + swig/python/install_python.cmake.in.
> Perhaps there's a better way...
>
> GDAL hopefully does the right thing for a typical distribution
> installation where CMAKE_INSTALL_PREFIX=/usr . For other prefixes, such as
> the implicit /usr/local prefix, you're on your own, and you will likely
> have to set the PYTHONPATH environment variable to point to
> /usr/local/lib/python3/dist-packages
>
> GDAL_PYTHON_INSTALL_PREFIX will not help you here (and I'm not sure to
> remember what the actual use case for it was)
>
> Even
> Le 08/03/2023 à 02:19, Michael Sumner a écrit :
>
> Hello, apologies as this is not so much a GDAL but a python question, but
> I'm interested in what others do, hopefully I'm missing a key step that's
> not hacky :)
>
> I made a gist to record the details:
>
> https://gist.github.com/mdsumner/526af876cfddaa5ff245ab376b3cec84
>
> The crux is,  GDAL has placed the python osgeo lib files into
> '/usr/local/lib/python3/dist-packages', but python is only looking in
> '/usr/local/lib/python3.10/dist-packages' (as seen in sys.path).
>
> Is there a config step with cmake, or a subsequent step with python config
> that I should be doing?
>
> Is it GDAL_PYTHON_INSTALL_PREFIX ?
>
> Thank you.
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] question about python installation

[Michael Sumner]

Hello, apologies as this is not so much a GDAL but a python question, but
I'm interested in what others do, hopefully I'm missing a key step that's
not hacky :)

I made a gist to record the details:

https://gist.github.com/mdsumner/526af876cfddaa5ff245ab376b3cec84

The crux is,  GDAL has placed the python osgeo lib files into
'/usr/local/lib/python3/dist-packages', but python is only looking in
'/usr/local/lib/python3.10/dist-packages' (as seen in sys.path).

Is there a config step with cmake, or a subsequent step with python config
that I should be doing?

Is it GDAL_PYTHON_INSTALL_PREFIX ?

Thank you.


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Convert NetCDF to GeoTIFF WGS 84

[Michael Sumner]

oh dear I missed a key element in the actual command:

gdalwarp band1.vrt  tiff_test.tif -ts 4000 4000 -te -116.000
660.000 4.000  780.000 -t_srs "EPSG:3857"

Best

On Thu, Feb 23, 2023 at 6:35 PM Michael Sumner  wrote:

> in short, take the grid specification you want and pass it to gdalwarp,
> you need the subdataset syntax of the particular variable
>
> you need target extent "-te xmin ymin xmax ymax" and  target size "-ts
> width height" (or use target resolution "-tr "  for a given pixel size)
> and -t_srs "EPSG:3857" (a compressed lookup syntax of the full srs string,
> see it at the bottom of the WKT)
>
> First check you have a fully fledged srs string in the source "Coordinate
> System is:"
>
> gdalinfo NETCDF:"/home/me/netcdf_file.nc":hs
>
> if not you can add it with -a_srs in the gdal_translate step below.
>
> Always a good idea to do a quick test first (take away the -b argument
> once you're satisfied to write all bands to the final tif).
>
> gdal_translate NETCDF:"/home/me/netcdf_file.nc":hs -b 1 band1.vrt
> gdalwarp band1.vrt  tiff_test.tif -ts 4000 4000 -te -116.000
> 660.000 4.000  780.000 -t_srs
>
> We can't advise on the right -a_srs if the "Coordinate System is:" element
> is missing in the source gdalinfo, please share that output if you have
> follow up questions.  For deeper investigation having the actual file can
> be necessary.
>
> Consider the -r argument for resampling by gdalwarp, nearest neighbour is
> the default.
>
> Finally, you say "WGS84" but provide Pseudo Mercator world grid -
> sometimes "WGS84" means "latitude/longitude ..." - but I've just taken your
> request as given. The principle always applies, good extent and srs info
> for the input, provide any extent, dimension, srs request for the output
> and the warper sorts it out.
>
> HTH, Mike
>
>
>
> On Wed, Feb 15, 2023 at 5:12 AM Javier Garcia 
> wrote:
>
>> Hi,
>>
>> I would like to convert the next NetCDF file to a GeoTIFF WGS 84.
>>
>> NetCDF file `gdalinfo` output:
>>
>> Driver: netCDF/Network Common Data Format
>> Files: /home/me/netcdf_file.nc
>> Size is 512, 512
>> Metadata:
>>   NC_GLOBAL#altimetry=altiCERSAT
>>   NC_GLOBAL#altitude_resolution=n/a
>>   NC_GLOBAL#area=Ireland 2 min wave grid
>>   NC_GLOBAL#BETAMAX=1.52
>>   NC_GLOBAL#CICE0=0.25
>>   NC_GLOBAL#CICEN=0.75
>>   NC_GLOBAL#comment=
>>   NC_GLOBAL#contact=wa...@mycompany.com
>>   NC_GLOBAL#Conventions=CF-1.5
>>   NC_GLOBAL#creation_date=2017-07-27T12:20:30Z
>>   NC_GLOBAL#data_type=SWASH straight grid
>>   NC_GLOBAL#easternmost_longitude=-1.028
>>   NC_GLOBAL#easting=longitude
>>   NC_GLOBAL#field_type=hourly
>>   NC_GLOBAL#FLAGTR=4
>>   NC_GLOBAL#forcing_ice=iceCFSR
>>   NC_GLOBAL#forcing_wind=wndCFSR
>>   NC_GLOBAL#forecast_range=analysis
>>   NC_GLOBAL#forecast_type=hindcast
>>   NC_GLOBAL#grid=irluk-2m
>>   NC_GLOBAL#grid_projection=n/a
>>   NC_GLOBAL#history=2017-07-27T12:20:30Z : Creation
>>   NC_GLOBAL#institution=My Company
>>   NC_GLOBAL#institution_references=http://www.mycompany.com/
>>   NC_GLOBAL#latitude_resolution=   0.033
>>   NC_GLOBAL#longitude_resolution=   0.033
>>   NC_GLOBAL#maximum_altitude=9000 m
>>   NC_GLOBAL#minimum_altitude=-12000 m
>>   NC_GLOBAL#netcdf_version=netCDF 4.3.0
>>   NC_GLOBAL#northernmost_latitude=59.500
>>   NC_GLOBAL#northing=latitude
>>   NC_GLOBAL#product_name=ww3.irluk-2m.199501_hs.nc
>>   NC_GLOBAL#product_version=1.0
>>   NC_GLOBAL#REFCOAST=0.05
>>   NC_GLOBAL#references=http://www.mycompany.com/waves/
>>   NC_GLOBAL#REFICEBERG=0.3
>>   NC_GLOBAL#REFSUBGRID=0.05
>>   NC_GLOBAL#run_time=2017-07-27T12:20:30Z
>>   NC_GLOBAL#SDS4 namelist parameter WHITECAPWIDTH=0.3001
>>   NC_GLOBAL#source=TWM WAVEWATCH III (R) v4.18
>>   NC_GLOBAL#southernmost_latitude=49.500
>>   NC_GLOBAL#start_date=1995-01-01T00:00:00Z
>>   NC_GLOBAL#stop_date=1995-01-31T23:00:00Z
>>   NC_GLOBAL#title=TWM SWASH HINDCAST
>>   NC_GLOBAL#WAVEWATCH_III_switches=F90 NOGRB NOPA LRB4 SHRD PR3 UQ FLX0
>> LN1 ST4 STAB0 NL1 BT4 DB1 MLIM TR0 BS0 IC0 REF1 XX0 WNT1 WNX1 CRT1 CRX1 O0
>> O1 O2 O3 O4 O5 O6 O7 011 NC4 TRKNC IS0
>>   NC_GLOBAL#WAVEWATCH_III_version_number=4.18b
>>   NC_GLOBAL#westernmost_longitude=-12.000
>> Subdatasets:
>>   SUBDATASET_1_NAME=NETCDF:"/home/me/netcdf_file.nc":MAPSTA
>>   SUBDATASET_1_DESC=[301x330] status map (16-bit integer)
>>

Re: [gdal-dev] Convert NetCDF to GeoTIFF WGS 84

[Michael Sumner]

t; DATUM["World Geodetic System 1984",
> ELLIPSOID["WGS 84",6378137,298.257223563,
> LENGTHUNIT["metre",1]]],
> PRIMEM["Greenwich",0,
> ANGLEUNIT["degree",0.0174532925199433]],
> ID["EPSG",4326]],
> CONVERSION["Popular Visualisation Pseudo-Mercator",
> METHOD["Popular Visualisation Pseudo Mercator",
> ID["EPSG",1024]],
> PARAMETER["Latitude of natural origin",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8801]],
> PARAMETER["Longitude of natural origin",0,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8802]],
> PARAMETER["False easting",0,
> LENGTHUNIT["metre",1],
> ID["EPSG",8806]],
> PARAMETER["False northing",0,
> LENGTHUNIT["metre",1],
> ID["EPSG",8807]]],
> CS[Cartesian,2],
> AXIS["easting (X)",east,
> ORDER[1],
> LENGTHUNIT["metre",1]],
> AXIS["northing (Y)",north,
> ORDER[2],
> LENGTHUNIT["metre",1]],
> USAGE[
> SCOPE["unknown"],
> AREA["World - 85°S to 85°N"],
> BBOX[-85.06,-180,85.06,180]],
> ID["EPSG",3857]]
> Data axis to CRS axis mapping: 1,2
> Origin = (-116.000,779.8137355)
> Pixel Size = (300.000,-300.000)
> Metadata:
>   AREA_OR_POINT=Area
> Image Structure Metadata:
>   COMPRESSION=LZW
>   INTERLEAVE=PIXEL
> Corner Coordinates:
> Upper Left  (-116.000, 780.000) ( 10d25'13.65"W, 57d11'40.60"N)
> Lower Left  (-116.000, 660.000) ( 10d25'13.65"W, 50d52'46.07"N)
> Upper Right (   4.000, 780.000) (  0d21'33.57"E, 57d11'40.60"N)
> Lower Right (   4.000, 660.000) (  0d21'33.57"E, 50d52'46.07"N)
> Center  ( -56.000, 720.000) (  5d 1'50.04"W, 54d 9'26.21"N)
> Band 1 Block=4000x1 Type=Byte, ColorInterp=Red
>   Mask Flags: PER_DATASET ALPHA
> Band 2 Block=4000x1 Type=Byte, ColorInterp=Green
>   Mask Flags: PER_DATASET ALPHA
> Band 3 Block=4000x1 Type=Byte, ColorInterp=Blue
>   Mask Flags: PER_DATASET ALPHA
> Band 4 Block=4000x1 Type=Byte, ColorInterp=Alpha
>
>
> Any help will be appreciated.
>
> Javier
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] help with (docker) build workflow

[Michael Sumner]

Ah, awesome - thanks Even and Daniel, both very helpful.

Cheers, Mike



On Mon, Feb 20, 2023 at 5:36 AM Even Rouault 
wrote:

> Michael,
>
> for your use case, you are only interested in the "builder" image, not the
> "runner" one.
>
> So just build it with something like: docker build . --target builder -t
> whatever_tag_name_you_want
>
> Even
> Le 19/02/2023 à 19:08, Michael Sumner a écrit :
>
> I would like to build GDAL in an image on the basis of a dockerfile used
> for CI, I'm confused by the layer process that copies the build outputs for
> the final steps
>
>
> https://github.com/OSGeo/gdal/blob/master/docker/ubuntu-small/Dockerfile#L231
>
> I want to be able to run that image and then do my own builds with cmake -
> but the layer COPY process doesn't include cmake or the other prereqs used
> above (to save space in the image).
>
> Is it easy to include the build capability as a layer-copy step?
>
> Or is this perhaps a bad idea, I'm looking for the easiest way to get a
> build system running with docker - I can unpick the builder/COPY process,
> but maybe I'm missing something easy.
>
> Thank you
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] help with (docker) build workflow

[Michael Sumner]

I would like to build GDAL in an image on the basis of a dockerfile used
for CI, I'm confused by the layer process that copies the build outputs for
the final steps

https://github.com/OSGeo/gdal/blob/master/docker/ubuntu-small/Dockerfile#L231

I want to be able to run that image and then do my own builds with cmake -
but the layer COPY process doesn't include cmake or the other prereqs used
above (to save space in the image).

Is it easy to include the build capability as a layer-copy step?

Or is this perhaps a bad idea, I'm looking for the easiest way to get a
build system running with docker - I can unpick the builder/COPY process,
but maybe I'm missing something easy.

Thank you

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] why padfExtent/s in gdal_alg.h ?

[Michael Sumner]

Thank you, gosh that does help explain a lot.

Cheers, Mike

On Wed, 8 Feb 2023, 21:42 Thomas Knudsen,  wrote:

> >  Ah ok, mostly I was just confused why it's not an error for the compiler
>
> The compiler only needs to know the type of the argument - in a header
> file, the identifier name following
> the type is for human information only
>
>>
>>
> Den ons. 8. feb. 2023 kl. 00.49 skrev Michael Sumner :
>
>>
>>
>> On Wed, Feb 8, 2023 at 10:44 AM Even Rouault 
>> wrote:
>>
>>> Michael,
>>>
>>> Feel free to submit a fix as there's obviously no justification for the
>>> discrepancy other than just being accidental.
>>>
>>
>> Ah ok, mostly I was just confused why it's not an error for the compiler
>> ... but I'm getting used to surprises learning C++ :)
>>
>> PR incoming, thanks!
>>
>> Cheers, Mike
>>
>>
>>
>>> There are likely hundreds of similar occurrences in the code base.  I
>>> believe there's a compiler warning (or at least cppcheck) to detect this,
>>> but it is not activated.
>>>
>>> Even
>>> Le 08/02/2023 à 00:35, Michael Sumner a écrit :
>>>
>>> This is not causing any problems afaict, but I happened to notice the
>>> plural for 'padfExtents' in the gdal_alg.h signature for
>>> GDALSuggestedWarpOutput2 compared to gdaltransformer.cpp
>>>
>>>
>>> https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdal_alg.h#L264
>>>
>>>
>>> https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdaltransformer.cpp#L412
>>>
>>>
>>> The discrepancy is visible in the doc here:
>>>
>>>
>>> https://gdal.org/api/gdal_alg.html#_CPPv424GDALSuggestedWarpOutput212GDALDatasetH19GDALTransformerFuncPvPdPiPiPdi
>>>
>>> Cheers, Mike
>>>
>>>
>>>
>>> --
>>> Michael Sumner
>>> Software and Database Engineer
>>> Australian Antarctic Division
>>> Hobart, Australia
>>> e-mail: mdsum...@gmail.com
>>>
>>> ___
>>> gdal-dev mailing 
>>> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>
>>> -- http://www.spatialys.com
>>> My software is free, but my time generally not.
>>>
>>>
>>
>> --
>> Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> Hobart, Australia
>> e-mail: mdsum...@gmail.com
>> ___
>> gdal-dev mailing list
>> gdal-dev@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] why padfExtent/s in gdal_alg.h ?

[Michael Sumner]

On Wed, Feb 8, 2023 at 10:44 AM Even Rouault 
wrote:

> Michael,
>
> Feel free to submit a fix as there's obviously no justification for the
> discrepancy other than just being accidental.
>

Ah ok, mostly I was just confused why it's not an error for the compiler
... but I'm getting used to surprises learning C++ :)

PR incoming, thanks!

Cheers, Mike



> There are likely hundreds of similar occurrences in the code base.  I
> believe there's a compiler warning (or at least cppcheck) to detect this,
> but it is not activated.
>
> Even
> Le 08/02/2023 à 00:35, Michael Sumner a écrit :
>
> This is not causing any problems afaict, but I happened to notice the
> plural for 'padfExtents' in the gdal_alg.h signature for
> GDALSuggestedWarpOutput2 compared to gdaltransformer.cpp
>
>
> https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdal_alg.h#L264
>
>
> https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdaltransformer.cpp#L412
>
>
> The discrepancy is visible in the doc here:
>
>
> https://gdal.org/api/gdal_alg.html#_CPPv424GDALSuggestedWarpOutput212GDALDatasetH19GDALTransformerFuncPvPdPiPiPdi
>
> Cheers, Mike
>
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] why padfExtent/s in gdal_alg.h ?

[Michael Sumner]

This is not causing any problems afaict, but I happened to notice the
plural for 'padfExtents' in the gdal_alg.h signature for
GDALSuggestedWarpOutput2 compared to gdaltransformer.cpp

https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdal_alg.h#L264

https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdaltransformer.cpp#L412


The discrepancy is visible in the doc here:

https://gdal.org/api/gdal_alg.html#_CPPv424GDALSuggestedWarpOutput212GDALDatasetH19GDALTransformerFuncPvPdPiPiPdi

Cheers, Mike



-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] clip raster to one cell at coordinates

[Michael Sumner]

Hi Karsten,

An option is to run gdallocationinfo, get the col (P) and row (L) value and
then feed that into

gdal_translate -srcwin P L col -outsize 1 1  in.tif out.tif ##I think
that's right

You otherwise have to do the "grid logic", which (while easy ...) I'd
rather like to see exposed explicitly in GDAL as functionality and I'm
exploring this in various ways atm.  (I think you can't use -projwin as it
seems to need at least a half pixel width and height to count as non zero
so it's still a bit tricky to get the right values).

Cheers, Mike


On Fri, Jan 27, 2023 at 5:01 PM karsten  wrote:

> Hi All,
>
> in an application I have a user click on a map and I need to
> programmatically clip a raster layer
> to the one raster cell for those coordinates at the click location.
> I was looking e.g. at https://gdal.org/programs/gdal_rasterize.html to
> run as a command
> line tool, but so far could not figure out if I can some can "feed" that a
> vector point to the tool as an input
> without first writing a shape file with the point onto disk.
>
> Does somebody have any pointers to a command line tool (or a python
> example) I could use to achieve this ?
>
> Thanks
> Karsten
>
> Karsten Vennemann
> Principal
>
> Terra GIS LTD
> 2119 Boyer Ave E
> Seattle, WA  98112
> USA
> www.terragis.net
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Call for review and comments on RFC90: Direct access to compressed raster data

[Michael Sumner]

beautiful!

This has been low level on my mind to ask about. Thanks!

Cheers, Mike

On Fri, 6 Jan 2023, 01:42 Even Rouault,  wrote:

> Hi,
>
> I've submitted RFC90: Direct access to compressed raster data for review
> and comments:
>
>  https://github.com/OSGeo/gdal/pull/7020
>
> Summary: """
>
> The document proposes 2 new methods to directly obtain the content of a
> window of interest of a raster dataset in its native compressed format.
> Those methods can be optionally implemented and used by drivers to perform:
> - extraction of a compressed tile as a standalone file from a container
> format (GeoTIFF, GeoPackage, etc.)
> - creation of mosaics from a set of tiles in individual files
> - lossless conversion between container formats using the same
> compression method
> """
>
> 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
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalwarp vs gdal_translate for scaling

[Michael Sumner]

ugh, sorry just realized my comment (sentence begins 'Further ...') about
"scaling" is completely irrelevant - you mean spatial and not data scaling.
Just ignore that one sentence.

Cheers, Mike



On Fri, Jan 6, 2023 at 10:36 AM Michael Sumner  wrote:

> in short, warp will give you exactly the target extent of a (possibly) new
> grid specification in that -te subdivided by the -ts.
>
> translate will give you approximately the extent provided to projwin, but
> snapped to the source grid alignment at a particular zoom implied by the
> projwin subdivided by the outsize.
>
> Further, translate can do scaling but warp cannot, which I guess is why
> you have a vrt in the workflow - you can scale a source virtually with
> translate to vrt, then warp (or translate) that.
>
> (Finally it's a bit sloppy to have such arbitrary extent values that warp
> will preserve and I routinely snap them to the source alignment, and I'm
> interested to provide such a feature in the warper options too.)
>
> Hope that helps!
>
>
>
>
>
>
> On Fri, 6 Jan 2023, 03:16 Daniel Mannarino, 
> wrote:
>
>> Hi folks!
>>
>>
>>
>> We have a process that generates web mercator zoom levels by resampling
>> from the next higher up, from the original resolution down to zoom level 0.
>> Currently we do so with gdalwarp, but a team member realized that
>> gdal_translate is capable of the same, and at first glance is much faster.
>> Unfortunately, using the two commands (below) which I intended to do
>> exactly the same thing winds up with slightly different results (can/should
>> I post a small image here to illustrate?).
>>
>> I won't ask "Which is correct?" because the answer is almost certainly
>> "That depends on what you want". But can someone help me understand the
>> difference between what is done with the gdalwarp command versus the
>> gdal_translate one so that I can decide if using gdal_translate is an
>> acceptable solution? Here are the two commands:
>>
>>
>>
>> gdalwarp -te 12523442.714243278 0.0 15028131.257091936 2504688.5428798534
>> -ts 65536 65536 -r mode -co TILED=YES ../source_tiles/everything.vrt
>> 007R_013C_warp_mode.tif
>>
>>
>>
>> gdal_translate -r mode -projwin 12523442.714243278 2504688.5428798534
>> 15028131.257091936 0.0 -outsize 65536 65536 -co TILED=YES
>> ../source_tiles/everything.vrt 007R_013C_translate_mode.tif
>>
>>
>>
>> Thank you!
>>
>> Daniel Mannarino
>> ___
>> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalwarp vs gdal_translate for scaling

[Michael Sumner]

in short, warp will give you exactly the target extent of a (possibly) new
grid specification in that -te subdivided by the -ts.

translate will give you approximately the extent provided to projwin, but
snapped to the source grid alignment at a particular zoom implied by the
projwin subdivided by the outsize.

Further, translate can do scaling but warp cannot, which I guess is why you
have a vrt in the workflow - you can scale a source virtually with
translate to vrt, then warp (or translate) that.

(Finally it's a bit sloppy to have such arbitrary extent values that warp
will preserve and I routinely snap them to the source alignment, and I'm
interested to provide such a feature in the warper options too.)

Hope that helps!






On Fri, 6 Jan 2023, 03:16 Daniel Mannarino, 
wrote:

> Hi folks!
>
>
>
> We have a process that generates web mercator zoom levels by resampling
> from the next higher up, from the original resolution down to zoom level 0.
> Currently we do so with gdalwarp, but a team member realized that
> gdal_translate is capable of the same, and at first glance is much faster.
> Unfortunately, using the two commands (below) which I intended to do
> exactly the same thing winds up with slightly different results (can/should
> I post a small image here to illustrate?).
>
> I won't ask "Which is correct?" because the answer is almost certainly
> "That depends on what you want". But can someone help me understand the
> difference between what is done with the gdalwarp command versus the
> gdal_translate one so that I can decide if using gdal_translate is an
> acceptable solution? Here are the two commands:
>
>
>
> gdalwarp -te 12523442.714243278 0.0 15028131.257091936 2504688.5428798534
> -ts 65536 65536 -r mode -co TILED=YES ../source_tiles/everything.vrt
> 007R_013C_warp_mode.tif
>
>
>
> gdal_translate -r mode -projwin 12523442.714243278 2504688.5428798534
> 15028131.257091936 0.0 -outsize 65536 65536 -co TILED=YES
> ../source_tiles/everything.vrt 007R_013C_translate_mode.tif
>
>
>
> Thank you!
>
> Daniel Mannarino
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Specify Grid Rotation for gdal.Warp() with Geolocation Arrays

[Michael Sumner]

yes I see, well you can probably set a target GeoTIFF with the geotransform
that you want - and gdalwarp to that as the target.  I don't believe you
can do that with gdalwarp directly, and I have no idea with python - and
finally I'm sorry to add irrelevant discussion - but this is interesting
and I'm going to try it out in R and C++ and would love to see a python
example too.

Cheers, Mike





On Sat, Nov 5, 2022 at 2:21 PM Brendan Heberlein 
wrote:

> Hi Mike,
>
> regardless of the map projection, it may be desirable to have the grid
> oriented in one direction or another, e.g. due to the spatial distribution
> of data points. Although the CRS/SRS designates the coordinate axes, these
> are not always aligned with the desired grid orientation.
>
> So, basically I want to be able to specify an arbitrary geotransform array
> for the output grid. This would allow the grid to be assigned an arbitrary
> orientation within the CRS frame.
>
> Thanks,
> Brendan
> ------
> *From:* Michael Sumner 
> *Sent:* Friday, November 4, 2022 6:57 PM
> *To:* Brendan Heberlein 
> *Cc:* gdal-dev 
> *Subject:* Re: [gdal-dev] Specify Grid Rotation for gdal.Warp() with
> Geolocation Arrays
>
> well, the map projection does this - there are many possibilities. Or, do
> you have another geolocation array you want as the target?
>
> If you have an example I'm happy to try a few things, but I'd assumed you
> had a target map projection in mind. What kind of crs to choose depends on
> your goal, or is purely to not have north up?
>
> Perhaps there's something terminology wise I'm missing ...
>
> Mike
>
> On Sat, 5 Nov 2022, 10:50 Brendan Heberlein,  wrote:
>
> Hi Mike,
>
> thanks for the response.
>
> Can you clarify how I would go about specifying the grid orientation for
> the output raster? None of -te, -ts or -t_srs address this.
>
>
> --
> *From:* Michael Sumner 
> *Sent:* Friday, November 4, 2022 6:41 PM
> *To:* Brendan Heberlein 
> *Cc:* gdal-dev 
> *Subject:* Re: [gdal-dev] Specify Grid Rotation for gdal.Warp() with
> Geolocation Arrays
>
> that's exactly what the warper does with geolocation arrays, set the
> target extent, dimension, and crs with -te, -ts, -t_srs with gdalwarp.
>
> python will have those analogous controls for the warper.
>
> Cheers, Mike
>
>
>
> On Sat, 5 Nov 2022, 01:23 Brendan Heberlein via gdal-dev, <
> gdal-dev@lists.osgeo.org> wrote:
>
> Hello,
>
> I would like to be able to warp an image dataset using a geolocation
> array, while specifying the grid to which the output dataset is sampled.
> Specifically, I want to be able to warp the dataset to a grid which is not
> oriented North-up.
>
> Can GDAL support this currently? Or, what is the likelihood this could be
> supported in the future?
>
> I primarily rely on the Python bindings, so a solution within that
> framework would be ideal for me.
>
> Thanks!
>  — Brendan
>
>
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Specify Grid Rotation for gdal.Warp() with Geolocation Arrays

[Michael Sumner]

well, the map projection does this - there are many possibilities. Or, do
you have another geolocation array you want as the target?

If you have an example I'm happy to try a few things, but I'd assumed you
had a target map projection in mind. What kind of crs to choose depends on
your goal, or is purely to not have north up?

Perhaps there's something terminology wise I'm missing ...

Mike

On Sat, 5 Nov 2022, 10:50 Brendan Heberlein,  wrote:

> Hi Mike,
>
> thanks for the response.
>
> Can you clarify how I would go about specifying the grid orientation for
> the output raster? None of -te, -ts or -t_srs address this.
>
>
> ----------
> *From:* Michael Sumner 
> *Sent:* Friday, November 4, 2022 6:41 PM
> *To:* Brendan Heberlein 
> *Cc:* gdal-dev 
> *Subject:* Re: [gdal-dev] Specify Grid Rotation for gdal.Warp() with
> Geolocation Arrays
>
> that's exactly what the warper does with geolocation arrays, set the
> target extent, dimension, and crs with -te, -ts, -t_srs with gdalwarp.
>
> python will have those analogous controls for the warper.
>
> Cheers, Mike
>
>
>
> On Sat, 5 Nov 2022, 01:23 Brendan Heberlein via gdal-dev, <
> gdal-dev@lists.osgeo.org> wrote:
>
> Hello,
>
> I would like to be able to warp an image dataset using a geolocation
> array, while specifying the grid to which the output dataset is sampled.
> Specifically, I want to be able to warp the dataset to a grid which is not
> oriented North-up.
>
> Can GDAL support this currently? Or, what is the likelihood this could be
> supported in the future?
>
> I primarily rely on the Python bindings, so a solution within that
> framework would be ideal for me.
>
> Thanks!
>  — Brendan
>
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Specify Grid Rotation for gdal.Warp() with Geolocation Arrays

[Michael Sumner]

that's exactly what the warper does with geolocation arrays, set the target
extent, dimension, and crs with -te, -ts, -t_srs with gdalwarp.

python will have those analogous controls for the warper.

Cheers, Mike



On Sat, 5 Nov 2022, 01:23 Brendan Heberlein via gdal-dev, <
gdal-dev@lists.osgeo.org> wrote:

> Hello,
>
> I would like to be able to warp an image dataset using a geolocation
> array, while specifying the grid to which the output dataset is sampled.
> Specifically, I want to be able to warp the dataset to a grid which is not
> oriented North-up.
>
> Can GDAL support this currently? Or, what is the likelihood this could be
> supported in the future?
>
> I primarily rely on the Python bindings, so a solution within that
> framework would be ideal for me.
>
> Thanks!
>  — Brendan
>
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Zarr syntax question

[Michael Sumner]

Thanks Even!

I THOUGHT I TRIED THAT  doh. Very much appreciated.

Cheers, Mike


On Thu, 8 Sept 2022, 19:39 Even Rouault,  wrote:

> Michael,
>
> This is both a usage issue, and the driver not erroring out on that
> invalid use whereas it should.
>
> You should prefix the URL with /vsicurl/ , and if you use Bash as a shell,
> to make sure the double quotes are preserved up to GDAL, you need to
> surround the whole string with single quote characters.
>
> So:
>
> gdalmdiminfo 'ZARR:"
> /vsicurl/https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr
> <https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr>
> "'
>
> and similarly for gdalinfo:
>
> gdalinfo  'ZARR:"/vsicurl/
> https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr;'
>
> gdalinfo 'ZARR:"/vsicurl/
> https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr":/precip:0'
> (0 here for the 2D slice at time index = 0)
>
> In https://github.com/OSGeo/gdal/pull/6330, I've improved doc and error
> messages.
>
> Even
>
>
> Le 08/09/2022 à 05:53, Michael Sumner a écrit :
>
> Hello, I'm trying to connect to this zarr:
>
> https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr
>
> I can get this sort of hull description with mdim, but only "cannot find
> group" with gdalinfo.
>
> gdalmdiminfo ZARR:"
> https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr;
> {
>   "type": "group",
>   "driver": "Zarr",
>   "name": "/"
> }
>
>
> gdalinfo ZARR:"
> https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr;
> ERROR 1: Cannot find group ncsa.osn.xsede.org
>
> That group error makes me think there's a kind of nested dataset issue to
> navigate?  But I can't see how zarr or gdal would expose those?   There is
> a variable called 'precip'.
>
> I'm comparing to the python in this notebook, which works fine on the same
> system. I'm using GDAL built from source.
>
>
> https://notebooksharing.space/view/c6c1f3a7d0c260724115eaa2bf78f3738b275f7f633c1558639e7bbd75b31456#displayOptions=\
>
> Cheers, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 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

[gdal-dev] Zarr syntax question

[Michael Sumner]

Hello, I'm trying to connect to this zarr:

https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr

I can get this sort of hull description with mdim, but only "cannot find
group" with gdalinfo.

gdalmdiminfo ZARR:"
https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr;
{
  "type": "group",
  "driver": "Zarr",
  "name": "/"
}


gdalinfo ZARR:"
https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr;
ERROR 1: Cannot find group ncsa.osn.xsede.org

That group error makes me think there's a kind of nested dataset issue to
navigate?  But I can't see how zarr or gdal would expose those?   There is
a variable called 'precip'.

I'm comparing to the python in this notebook, which works fine on the same
system. I'm using GDAL built from source.

https://notebooksharing.space/view/c6c1f3a7d0c260724115eaa2bf78f3738b275f7f633c1558639e7bbd75b31456#displayOptions=\

Cheers, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] EXT: Re: Rotated lat lon projection

[Michael Sumner]

Ok, thanks for the file. Here's a gist just to show you it works, I use R
so I've got some checks I do around the core gdalwarp call which can be as
minimal as

gdalwarp  20220902T00Z_MSC_HRDPS_HGT_SFC_RLatLon0.0225_P001.grib2 out.tif
-t_srs "OGC:CRS84"

You should consider the specification of the target grid yourself, with a
combination of some of -tr, -ts, -te, and -t_srs arguments. Also the -r
resampling argument, default is nearest neighbour.

https://gist.github.com/mdsumner/d73748c44aef4a0923760401d3dc605d

Cheers, Mike


On Mon, Sep 5, 2022 at 5:21 PM De Juan De España, Javier <
jajuan_...@acciona.com> wrote:

> Hi Mike!
>
>
>
> Thank you very much for your response.
>
>
>
> Attached you can find a grib2. Trying to test this on QGIS (Uses GDAL
> behind), it automatically opens it with  right projection so I guess GDAL
> Will do the same.
>
>
>
> In this file, first grid point should be something like:
>
>
>
> Thank you very much!!
>
>
>
>
>
> *De:* Michael Sumner 
> *Enviado el:* sábado, 3 de septiembre de 2022 8:59
> *Para:* De Juan De España, Javier 
> *CC:* gdal-dev 
> *Asunto:* EXT: Re: [gdal-dev] Rotated lat lon projection
>
>
>
> First step is try shoving it through the warper:
>
>
>
> gdalwarp in.nc
> <https://clicktime.symantec.com/15t5ZrhtFMk6q5dh8yXpr?h=f3NsNfYLCjuppp90Nv9-oEa-gZU0oRBfBuHIbBwMVFs==http://in.nc>
> out.tif
>
>
>
> set lonlat extent with -te and resolution with -ts or -tr
>
>
>
> GDAL knows about projections and internal coordinate arrays and often this
> is enough, the smart things come down to what output grid you define and
> what resampling algorithm is used and whether it's appropriate (beware of
> vector quantities for example).
>
>
>
> Use a VRT intermediate to first select out one band to test with
> (gdal_translate in.nc
> <https://clicktime.symantec.com/15t5ZrhtFMk6q5dh8yXpr?h=f3NsNfYLCjuppp90Nv9-oEa-gZU0oRBfBuHIbBwMVFs==http://in.nc>
> in1.vrt -b 1). If you have a test file I'm happy to explore.
>
>
>
> Even if the auto stuff doesn't work you can engineer a lot with just
> command line args and VRT even for complex grids.
>
>
>
> Cheers, Mike
>
>
>
> On Thu, 1 Sept 2022, 23:57 De Juan De España, Javier, <
> jajuan_...@acciona.com> wrote:
>
> Hi!
>
>
>
> We are developping a meteorology project in the Canadian domain where we
> use gdal to manage grib files.
>
> These *files are “rotated lat-lon” and we would like to switch to
> “regular lat-lon”* easily. I’m concerned that lastest versions of gdal
> can do this, but i would really appreciate some indications or information
> about how to do this.
>
>
>
> Thank very much!
>
>
>
> Javier de Juan
>
>
>
>
>
> [image: logo-bau]
>
>
>
>
>
>
>
> *Javier de Juan de España*
>
> Project Manager
>
> Smart Society Skill Center
>
> Digital Hub
>
> Avda Europa, 18
>
> 28108 Alcobendas (Madrid) – Spain
>
> +34 687 53 86 53
>
> [image: flecha]
>
> WWW.ACCIONA.COM
> <https://clicktime.symantec.com/15t5pMHjdCnt4v7TmejGi?h=bsqBMl9kpV7pNEC9p7lA_BzAvITG8lgP7QbVur0i_lg==https://www.acciona.com/es/>
>
> [image: twitter]
> <https://clicktime.symantec.com/15t5jX6TAb7HeyHYE6L86?h=-cIW40zXsiDGS34xBnF2qp0NFHQYcOLzae3_uhhM6P0==https://twitter.com/ACCIONA>
>
> [image: facebook]
> <https://clicktime.symantec.com/15t5uBV25pUUUrwPKD8RL?h=54T_VUKmdZlQ25O0tl3wf8KTqvyNg3vKMSd6fc0UwDU==https://www.facebook.com/acciona>
>
> [image: link]
> <https://clicktime.symantec.com/15t64qsb13qfJkbEQKvia?h=qP3wZRbx16RE2f9KqOZAFBfB0BsCiaiJdhQ6uIgX_GY==https://www.linkedin.com/company/acciona>
>
> [image: you]
> <https://clicktime.symantec.com/15t69g4sTfXFihR9wtKsC?h=smitgw2XvRmwLj9vHPFdbXjFyts7YW6yT7TAiJN1aOo==https://www.youtube.com/user/interacciona1?sub_confirmation%3D1>
>
> [image: insta]
> <https://clicktime.symantec.com/15t5z1gJYSA4tomJrmXZx?h=lsjgYLW3UgtJbbPBFBcqNBMP-9OUhZJgVLPeWqFyTaE==https://www.instagram.com/acciona/>
>
>
>
>
>
>
>
>
>
>
>
> Este mensaje se dirige exclusivamente a su destinatario, y puede contener
> información confidencial sometida a secreto profesional, o cuya divulgación
> esté legalmente prohibida. Cualquier opinión en él contenida es exclusiva
> de su autor y no representa necesariamente la opinión de la empresa. Si ha
> recibido este mensaje por error, le rogamos nos lo comunique de forma
> inmediata por esta misma vía y proceda a su eliminación, así como a la de
> cualquier documento adjunto al mismo. El correo electrónico vía Internet no
> es seguro y no se puede garantizar que no haya errores ya

Re: [gdal-dev] Rotated lat lon projection

[Michael Sumner]

First step is try shoving it through the warper:

gdalwarp in.nc out.tif

set lonlat extent with -te and resolution with -ts or -tr

GDAL knows about projections and internal coordinate arrays and often this
is enough, the smart things come down to what output grid you define and
what resampling algorithm is used and whether it's appropriate (beware of
vector quantities for example).

Use a VRT intermediate to first select out one band to test with
(gdal_translate in.nc in1.vrt -b 1). If you have a test file I'm happy to
explore.

Even if the auto stuff doesn't work you can engineer a lot with just
command line args and VRT even for complex grids.

Cheers, Mike

On Thu, 1 Sept 2022, 23:57 De Juan De España, Javier, <
jajuan_...@acciona.com> wrote:

> Hi!
>
>
>
> We are developping a meteorology project in the Canadian domain where we
> use gdal to manage grib files.
>
> These *files are “rotated lat-lon” and we would like to switch to
> “regular lat-lon”* easily. I’m concerned that lastest versions of gdal
> can do this, but i would really appreciate some indications or information
> about how to do this.
>
>
>
> Thank very much!
>
>
>
> Javier de Juan
>
>
>
>
>
> [image: logo-bau]
>
>
>
>
>
>
>
> *Javier de Juan de España*
>
> Project Manager
>
> Smart Society Skill Center
>
> Digital Hub
>
> Avda Europa, 18
>
> 28108 Alcobendas (Madrid) – Spain
>
> +34 687 53 86 53
>
> [image: flecha]
>
> WWW.ACCIONA.COM 
>
> [image: twitter] 
>
> [image: facebook] 
>
> [image: link] 
>
> [image: you]
> 
>
> [image: insta] 
>
>
>
>
>
>
>
>
>
>
>
> Este mensaje se dirige exclusivamente a su destinatario, y puede contener
> información confidencial sometida a secreto profesional, o cuya divulgación
> esté legalmente prohibida. Cualquier opinión en él contenida es exclusiva
> de su autor y no representa necesariamente la opinión de la empresa. Si ha
> recibido este mensaje por error, le rogamos nos lo comunique de forma
> inmediata por esta misma vía y proceda a su eliminación, así como a la de
> cualquier documento adjunto al mismo. El correo electrónico vía Internet no
> es seguro y no se puede garantizar que no haya errores ya que puede ser
> interceptado, modificado, perdido o destruido, o contener virus. Cualquier
> persona que se ponga en contacto con nosotros por correo electrónico se
> considerará que asume estos riesgos.
>
>
>
> This e-mail is addressed exclusively to the recipient and may contain
> privileged information under a professional confidential agreement or it
> may be against the law to disclose its contents. Any opinion contained in
> it belongs exclusively to his/her author and does not necessarily reflect
> the company’s view. If you receive this e-mail in error, please let us know
> immediately (by return e-mail) and proceed to its destruction, as well as
> any document attached to it. The sending of e-mails through the Internet is
> not safe and, therefore, error-free communications cannot be guaranteed, as
> they can be intercepted, changed, misled or destroyed or they might contain
> a virus. Any user contacting us through e-mails shall be understood to be
> assuming these risks.
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] question for scaling default

[Michael Sumner]

Hello, I have a new driver candidate, it's natively a type Byte and is a
percentage (or fraction) floating point value when unscaled.  I have a
couple of questions. .

https://github.com/OSGeo/gdal/pull/6183

I'd like to have it unscale by default, with band type Float32 - but be
able to opt out of scaling, and return the raw Byte type. I think that is a
good default since most usage will want the fraction value, but
occasionally useful to obtain the raw codes that include some special out
of range values (>255).

I think I can work out how to do this, but is there a driver that already
provides this kind of behaviour that I can learn from?

Also, I feel like fraction (0, 1) is a more sensible data value than
percentage (0, 100), although I've used percentage in the past with these
data. Any thoughts, or examples that exist in GDAL?

Thanks very much, Mike



-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] determine default SQL dialect

[Michael Sumner]

very helpful, thank you very much.  Quite a bit more to it than I thought.

Best, Mike



On Tue, Aug 2, 2022 at 7:27 PM Rahkonen Jukka <
jukka.rahko...@maanmittauslaitos.fi> wrote:

> Hi,
>
>
>
> About INDIRECT_SQLITE dialect, it is documented in
> https://gdal.org/user/sql_sqlite_dialect.html
>
> “If the datasource is SQLite database (GeoPackage, SpatiaLite) then SQLite
> dialect acts as native SQL dialect and Virtual Table Mechanism is not used.
> It is possible to force GDAL to use Virtual Tables even in this case by
> specifying “-dialect INDIRECT_SQLITE”. This should be used only when
> necessary, since going through the virtual table mechanism might affect
> performance.”
>
> There is also an usage example for one special use case.
>
>
>
> ES is documented in https://gdal.org/drivers/vector/elasticsearch.html.
>
> MONGO_DB is documented in https://gdal.org/drivers/vector/mongodbv3.html.
>
>
>
> -Jukka Rahkonen-
>
>
>
> *Lähettäjä:* gdal-dev  *Puolesta *Michael
> Sumner
> *Lähetetty:* tiistai 2. elokuuta 2022 10.43
> *Vastaanottaja:* gdal-dev 
> *Aihe:* [gdal-dev] determine default SQL dialect
>
>
>
> Is it possible to determine the default SQL dialect  that will be used in
> the C++ API?
>
>
>
> It seems like the logic to pivot on dialect choice is only available to
> drivers internally (they know if they have one), so a map of formats is
> enough - but I wonder if I'm missing an API hook.
>
>
>
> Also, there appears to be dialect values DEBUG, INDIRECT_SQLITE, MONGO_DB,
> and ES in use but I gather they are for internal use (and so properly don't
> appear in the documentation).
>
>
>
> Thanks very much, Mike
>
>
>
>
> --
>
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] determine default SQL dialect

[Michael Sumner]

Excellent, thanks very much.

Best, Mike

On Wed, 3 Aug 2022, 02:30 Even Rouault,  wrote:

>
> Le 02/08/2022 à 18:25, Even Rouault a écrit :
> >
> > Le 02/08/2022 à 09:42, Michael Sumner a écrit :
> >> Is it possible to determine the default SQL dialect  that will be
> >> used in the C++ API?
> >
> > Short answer: no
> >
> > Long answer: for a driver that is backed by a SQL engine (SQLite,
> > GPKG, PG, MySQL, OCI, etc.), it will default to the RDBMS SQL engine,
> > otherwise to OGR SQL.
>
> Potential not so pretty way of knowing if the default SQL dialect is OGR
> SQL:
>
> run ExecuteSQL( "SELECT OGR_GEOM_WKT FROM {some_layer_name}" )
>
> if that works, then at 99.99% , you have OGR SQL.
>
> 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

[gdal-dev] determine default SQL dialect

[Michael Sumner]

Is it possible to determine the default SQL dialect  that will be used in
the C++ API?

It seems like the logic to pivot on dialect choice is only available to
drivers internally (they know if they have one), so a map of formats is
enough - but I wonder if I'm missing an API hook.

Also, there appears to be dialect values DEBUG, INDIRECT_SQLITE, MONGO_DB,
and ES in use but I gather they are for internal use (and so properly don't
appear in the documentation).

Thanks very much, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] default for RelativeToVRT field

[Michael Sumner]

On Sat, Apr 23, 2022 at 7:37 PM Even Rouault 
wrote:

> Michael,
>
> > Or is it just a matter of forcing an input path to be absolute?
>
> for regular use cases (ie VRT name not set to empty), you might not even
> been able to do that in master since
>
> https://github.com/OSGeo/gdal/commit/6d7001656d06af138f8e46ff1651056626db11e8
> which tries to write relative filename as much as possible, even if all
> paths are absolute
>
> Testing with your use case with empty VRT filename, I see the behavior
> has not changed, and that you'll indeed get relativeToVRT="0" only if
> using a source dataset with an absolute path. For most use cases where
> you use an empty filename, it doesn't matter much because you don't
> actually re-open it and use it immediately, and thus the source dataset
> is not re-opened. I can imagine though that you could get into trouble
> if serializing it from the xml:VRT content and re-opening the resulting
> VRT.
>
> > Is it possible to force GDAL to make the SourceFilename absolute by
> > setting options?
>
> The solution would probably be to add an explicit creation option to the
> VRT driver, like SOURCE_PATH=ABSOLUTE, that would override the current
> heuristics.
>
>

Thanks a lot Even, very helpful as always. I can see why relative is the
(hard) default, (for auxiliary files to move around with sources).

I do think I would benefit from this creation option to make some workflows
easier, but I still need a bit longer exploring to get all the parts clear
(to me) -  I'll  PR this at some point.

Best, Mike



> Even
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Hints that a raster band represents a DEM?

[Michael Sumner]

read from a global dem via the warper and do a correlation test?

:)


On Sat, Apr 30, 2022 at 11:11 AM Nyall Dawson 
wrote:

> Hi list,
>
> I'm looking for any tips on potential approaches I could use to
> automagically guess that a particular band from a data source
> represents an elevation surface.
>
> I haven't been able to find any common metadata components in the
> files I've investigated which could help indicate this, but maybe I'm
> missing something. Right now the best approach I can think of would be
> to test if the layer's CRS has a vertical component (but unfortunately
> even among DEM layers I think this will be rare!).
>
> Any suggestions?
> Nyall
> ___
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] default for RelativeToVRT field

[Michael Sumner]

I would like to control the relativeToVRT field by setting it to 0/false,
it seems the only way to do that is to give a non-local or absolute path to
the *output* file, or give an absolute path for the input. .

intif=autotest/gdrivers/data/gtiff/byte_signed.tif
abtif="$(pwd)"/$intif
gdal_translate $intif -of VRT here.vrt
gdal_translate $intif -of VRT /tmp/there.vrt
gdal_translate $abtif -of VRT here_abs.vrt

grep relative here.vrt
  a ...
grep relative /tmp/there.vrt
  /p ...
grep relative here_abs.vrt
  /p ...

It's the same in C++, if 'pszFilename' is relative the value is 1, if
absolute it is 0.

  GDALDataset *poSrcDS = (GDALDataset *)GDALOpenShared(pszFilename,
GA_ReadOnly);

  GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
  GDALDataset *poVRTDS = poDriver->CreateCopy("", poSrcDS, false, NULL,
NULL, NULL);

  const char *xmlvrt = poVRTDS->GetMetadata("xml:VRT")[0];

Is it possible to force GDAL to make the SourceFilename absolute by setting
options? Or is it just a matter of forcing an input path to be absolute?

I'm concerned about cases of relative paths with prefixes, I want to create
in-memory VRT text that *always* has absolute paths for exactly this
CreateCopy("", ) situation with no actual VRT file.

Thank you, Mike



--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL warp tutorial question (C vs. C++)

[Michael Sumner]

Excellent - I'll submit a PR for the include in the first snippet.

I've become more au fait with the different idioms today, and this
description is very helpful.

Best, Mike

On Tue, Apr 5, 2022 at 8:47 PM Even Rouault 
wrote:
>
>
> Le 05/04/2022 à 01:28, Michael Sumner a écrit :
> > Hello, the warp tutorial  needs at least inclusion of "cpl_conv.h" for
> > CPLMalloc to work.
> >
> > https://gdal.org/tutorials/warp_tut.html
> >
> > But, I'd like to ask why it uses the style of the C examples (as per
> > the Raster API tutorial
> > https://gdal.org/tutorials/raster_api_tut.html) rather than the C++ -
> > i.e. it uses 'GDALDatasetH hSrcDS' declarations rather than
> > 'GDALDataset *poDataset', and all the 'GDAL' rather than
> > '->' pointer syntax.
> >
> > Is that recommended, or is it just old style in this tutorial (because
> > the apps were originally C)?
>
> Probably just the preference of the author of the example. At least the
> first code snippet
> (https://gdal.org/tutorials/warp_tut.html#a-simple-reprojection-case)
> uses some C++ with the GDALWarpOperation class so it could be
> potentially converted to full C++ (but you'll need to keep CPLMalloc()
> and not use new[] because GDALDestroyWarpOptions() expects panSrcBands
> and panDstBands to have been allocated with CPLMalloc())
>
> Actually I see the second snippet
> (https://gdal.org/tutorials/warp_tut.html#creating-the-output-file) does
> use some bits of C++ with the OGRSpatialReference class. It could be
> either fully converted to C++ (but some functions like
> GDALCreateGenImgProjTransformer() or GDALSuggestedWarpOutput() only
> exist in the C API, so you'd have to cast C++ dataset objects to C
> handles at some point. This also applies to the first snippet, which
> explains that using C functions is more practical to avoid those casts))
> or C.
>
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
>


--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL warp tutorial question (C vs. C++)

[Michael Sumner]

Hello, the warp tutorial  needs at least inclusion of "cpl_conv.h" for
CPLMalloc to work.

https://gdal.org/tutorials/warp_tut.html

But, I'd like to ask why it uses the style of the C examples (as per the
Raster API tutorial https://gdal.org/tutorials/raster_api_tut.html) rather
than the C++ - i.e. it uses 'GDALDatasetH hSrcDS' declarations rather than
'GDALDataset *poDataset', and all the 'GDAL' rather than '->'
pointer syntax.

Is that recommended, or is it just old style in this tutorial (because the
apps were originally C)?

(I still get rather stuck on some of these issues, my C++ is rudimentary -
certainly I'm keen to contribute to docs and more once I understand
better).

Thanks a lot.

Best, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] regular grid longlat netcdf without crs meta

[Michael Sumner]

Thanks Even, that's excellent.

I agree about the providers, and I inform them (the OISST project replied
with a promise to fix in future. I wish they and many others would catch on
to how terrible NetCDF is for 2D regular grids and just use good geotiffs
but hoo boo whoa).

Extending GDAL has more positive impact (for me at least) as the heuristics
will catch many other products.

Best, Mike


On Tue, 15 Mar 2022, 22:38 Even Rouault,  wrote:
>
> Michael,
>
> Le 15/03/2022 à 06:06, Michael Sumner a écrit :
>
> Hi,
>
> these files don't present to GDAL with a CRS, but it seems not out of
place to assume that they can use 'oSRS.SetWellKnownGeogCS("WGS84")', when
no other indication is found and there are dimensions 'lon,lat' or
'longitude,latitude'
>
>
https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/
>
> An example sds specifically:
>
> /vsicurl/NETCDF:"/vsicurl/
https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202202/oisst-avhrr-v02r01.20220201.nc
":sst
>
> Would that be reasonable? Perhaps controlled by an env var? Is there an
established control for this I'm missing? (I'm not looking for '-a_srs' or
VRT, I want the auto-detection to work for these and many other similar
file sources without augmentation).
>
> I wouldn't expose WGS84 by default, because it might be another Earth
geographic CRS, or even a non-Earth geographic CRS.
>
> Potential solution:
>
> - add a GDAL_NETCDF_DEFAULT_GEOGRAPHIC_CRS configuration option (that
could be set to EPSG:4326 e.g) that would be used in the situation where
the netCDF CF attributes to define a geographic CRS, but without the actual
CRS definition, are not set.
>
> - or, more involved but more powerful, add a configuration file in data/
, let's say a netcdf_product_conf.json or whatever, where you would define
a set of product types, with matching global attributes and their values
and the corresponding CRS.
>
> e.g (not sure my set of "matching_attributes" is appropriate for your use
case)
>
> {
> "product_type": {
> "noaa_ncei_oisst": {  // the id of the product here is completely
arbitrary
> "defining_global_attributes": {
> "title": "NOAA/NCEI 1/4 Degree Daily Optimum
Interpolation Sea Surface Temperature (OISST) Analysis, Version 2.1 -
Final",
> "geospatial_lat_units": "degrees_north",
> "geospatial_lon_units": "degrees_east"
> },
> "default_values": {
> "crs": "EPSG:4326"
> }
> }
> }
>
> }
>
> - or, both the easiest and hardest one, convince the data producer to
just use the netCDF CF conventions to their full capabilities.
>
> Even
>
>
> (It's not out of the question that I could submit a PR for this, just
flagging discussion in case I'm missing something.)
>
> Thanks.
>
> Best, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> --
> 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

Re: [gdal-dev] regular grid longlat netcdf without crs meta

[Michael Sumner]

gosh, very sorry I messed up the SDS string - I've posted as a gist here
with the gdalinfo output:

https://gist.github.com/mdsumner/615441365e3a0691c457493fe62b78bf

Best, Mike


On Tue, Mar 15, 2022 at 4:06 PM Michael Sumner  wrote:
>
> Hi,
>
> these files don't present to GDAL with a CRS, but it seems not out of
place to assume that they can use 'oSRS.SetWellKnownGeogCS("WGS84")', when
no other indication is found and there are dimensions 'lon,lat' or
'longitude,latitude'
>
>
https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/
>
> An example sds specifically:
>
> /vsicurl/NETCDF:"/vsicurl/
https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202202/oisst-avhrr-v02r01.20220201.nc
":sst
>
> Would that be reasonable? Perhaps controlled by an env var? Is there an
established control for this I'm missing? (I'm not looking for '-a_srs' or
VRT, I want the auto-detection to work for these and many other similar
file sources without augmentation).
>
> (It's not out of the question that I could submit a PR for this, just
flagging discussion in case I'm missing something.)
>
> Thanks.
>
> Best, Mike
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com



--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] regular grid longlat netcdf without crs meta

[Michael Sumner]

Hi,

these files don't present to GDAL with a CRS, but it seems not out of place
to assume that they can use 'oSRS.SetWellKnownGeogCS("WGS84")', when no
other indication is found and there are dimensions 'lon,lat' or
'longitude,latitude'

https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/

An example sds specifically:

/vsicurl/NETCDF:"/vsicurl/
https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202202/oisst-avhrr-v02r01.20220201.nc
":sst

Would that be reasonable? Perhaps controlled by an env var? Is there an
established control for this I'm missing? (I'm not looking for '-a_srs' or
VRT, I want the auto-detection to work for these and many other
similar file sources without augmentation).

(It's not out of the question that I could submit a PR for this, just
flagging discussion in case I'm missing something.)

Thanks.

Best, Mike


-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev