[gdal-dev] Solved: Re: Having a problem getting VRT ComplexSource to Scale

2020-01-23 Thread Stephen Woodbridge 
Never mind! I just checked the histogram and the data is getting scaled, 
I forgot that I was also scaling it in my mapfile so the results  where 
getting masked.


-Steve

On 1/23/2020 8:20 PM, Stephen Woodbridge wrote:

Hi all,

I'm having trouble getting my VRT file to scale the source data into 
my ColorTable.

I've been looking at this https://gdal.org/drivers/raster/vrt.html

The source GeoTiff has:

Band 1 Block=4500x1 Type=Float32, ColorInterp=Gray
  Minimum=1.082, Maximum=46.322, Mean=34.181, StdDev=2.049
  NoData Value=-3
  Overviews: 2250x2126, 1125x1063, 563x532, 282x266, 141x133, 71x67
  Metadata:
    STATISTICS_MAXIMUM=46.321998596191
    STATISTICS_MEAN=34.181107349018
    STATISTICS_MINIMUM=1.0819988250732
    STATISTICS_STDDEV=2.0492649377546

My goal is to scale the source 0-50 into 10-240 in the ColorTable.


The VRT looks like:


   -180.0 ,0.08 , 0 , 90.0 , 0 , -0.04 
  
  
  
    Palette
    
  
  
...
  
    
    
  relativeToVRT="1">HYCOM_tomorrow_salinity_0.tif

  1
  0
  1
  0
  
  
    
  


50 - 0 = 50
240 - 10 = 230
230 / 50 = 4.6 = ScaleRatio
10 = ScaleRatio

I've tried a lot of other value trying to play with the scaling.

I also tried using:

0

0
50
10
240 but that doesn't appear to work either. Using: 
GDAL 2.2.3, released 2017/11/20 on Unbuntu 18.04 My conclusion is that 
I'm not understanding the meaning of these fields or these are meant 
for something other than my goal. Any help would be appreciated. 
Thanks, -Steve




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

Re: [gdal-dev] Querying many raster pixels efficiently from Python

2020-01-23 Thread Patrick Young 
Putting something like mapserver/mapcache in front of your requests might
work, if I understand correctly.  We serve COGs out of s3 via WMS like this
and the performance is pretty nice.

See

https://github.com/pedros007/mapserver-docker

for some discussion on such a setup.

Best,
Patrick


On Thu, Jan 23, 2020, 6:44 AM Daniel Evans  wrote:

> Hi,
>
>
>
> I have a large (global, 30m resolution, 50GB+) GeoTIFF dataset, from which
> I need to read many (millions) of pixel values at given input coordinates.
> I’ve got reasonable performance out of the code, about a million queries
> over five minutes, but:
>
>
>
>1. There are actually twelve separate datasets of this size to query,
>not just one, so it takes approximately an hour.
>2. This is by far the slowest portion of the program, and the users
>demand speed!
>3. The users would also like to move towards higher resolution
>datasets, which we see run about 5x slower.
>4. When querying the data on a particular piece of network storage
>mounted as part of the local filesystem, we see a slowdown approaching two
>orders of magnitude – bulk file copies off the network storage are
>reasonable, but each IO request shows a significant overhead (up to a
>second), and GDAL is sending one for each coordinate queried.
>
>
>
> The implementation is in Python, directly calling down to GDAL. The short,
> long-running snippet of code which performs the actual queries the dataset,
> having converted real-world coordinates to pixels, is:
>
>
>
> value_arrays = (
>
> raster_ds.ReadAsArray(
>
> xoff=coord[0] - buffer_size,
>
> yoff=coord[1] - buffer_size,
>
> xsize=npix,
>
> ysize=npix
>
> ) for coord in offsets
>
> )
>
>
>
> There are a few things that are probably worth noting:
>
>1. It is not necessarily a single pixel that is being read – for each
>coordinate, the program may be asked to get all pixel values within a given
>radius (typically a couple of pixels), and use some function to summarise
>these into a single value (mean, median, …). GDAL currently returns a numpy
>array for each query, which is passed to the user-specified function after
>the snippet above.
>2. The dataset is made up of 2048x2048 LZW-Compressed tiles containing
>floating point data (essentially conforming to COG, but with no overviews),
>grouped together in a VRT (performance is identical with plain GeoTIFFs,
>though).
>
>
>
> Multiprocessing has not been found to help - we actually lose throughput
> as the disk read head is moving back and forth constantly. Better hardware
> (especially SSDs) is known to help, but no one wants to pay for that.
>
>
>
> We see no particular performance difference from setting
> GDAL_DISABLE_READDIR_ON_OPEN=TRUE, and GDAL_CACHEMAX is left at the default
> 5% (64GB+ RAM available).
>
>
>
>
>
>
>
> Does the Python interface to GDAL provide a way to supply a large number
> of offsets and get blocks of pixels back, avoiding the need to come back up
> to Python after each query? (I suspect not)
>
> Is there some way to optimise GDAL so that queries of files on the mounted
> network storage are more efficient?
>
>
>
>
>
>
>
> *Dr. Daniel Evans*
>
> Software Developer
>
>
>
> *Skype*
>
>
>
> *T* +44 (0) 1756 799919
> www.jbarisk.com
>
> [image: Visit our website]   [image: Follow us on
> Twitter] 
>
> Our postal address and registered office is JBA Risk Management Limited,
> 1 Broughton Park, Old Lane North, Broughton, Skipton, North Yorkshire BD23
> 3FD.
>
> *Find out more about us here: www.jbarisk.com 
> and **follow us on Twitter @JBARisk  and
> LinkedIn
> 
> *
>
> The JBA Group supports the JBA Trust.
>
> All JBA Risk Management's email messages contain confidential information
> and are intended only for the individual(s) named. If you are not the named
> addressee you should not disseminate, distribute or copy this e-mail.
> Please notify the sender immediately by email if you have received this
> email by mistake and delete this email from your system.
>
>
> JBA Risk Management Limited is registered in England, company number
> 07732946, 1 Broughton Park, Old Lane North, Broughton, Skipton, North
> Yorkshire, BD23 3FD, Telephone: +441756799919
>
>
> ___
> 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] Having a problem getting VRT ComplexSource to Scale

2020-01-23 Thread Stephen Woodbridge 

Hi all,

I'm having trouble getting my VRT file to scale the source data into my 
ColorTable.

I've been looking at this https://gdal.org/drivers/raster/vrt.html

The source GeoTiff has:

Band 1 Block=4500x1 Type=Float32, ColorInterp=Gray
  Minimum=1.082, Maximum=46.322, Mean=34.181, StdDev=2.049
  NoData Value=-3
  Overviews: 2250x2126, 1125x1063, 563x532, 282x266, 141x133, 71x67
  Metadata:
    STATISTICS_MAXIMUM=46.321998596191
    STATISTICS_MEAN=34.181107349018
    STATISTICS_MINIMUM=1.0819988250732
    STATISTICS_STDDEV=2.0492649377546

My goal is to scale the source 0-50 into 10-240 in the ColorTable.


The VRT looks like:


   -180.0 ,0.08 , 0 , 90.0 , 0 , -0.04 
  
  
  
    Palette
    
  
  
...
  
    
    
  relativeToVRT="1">HYCOM_tomorrow_salinity_0.tif

  1
  0
  1
  0
  
  
    
  


50 - 0 = 50
240 - 10 = 230
230 / 50 = 4.6 = ScaleRatio
10 = ScaleRatio

I've tried a lot of other value trying to play with the scaling.

I also tried using:

0

0
50
10
240 but that doesn't appear to work either. Using: GDAL 2.2.3, 
released 2017/11/20 on Unbuntu 18.04 My conclusion is that I'm not 
understanding the meaning of these fields or these are meant for 
something other than my goal. Any help would be appreciated. Thanks, -Steve


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

Re: [gdal-dev] OGR FileGDB driver: 'OBJECTID' not recognised as an available field

2020-01-23 Thread Tamas Szekeres 
Hi,

I've already reported an issue a while back related to the OGRSQL dialect
along with the FDGB driver, when the where clause contains expression for
the FID column. I'm still experiencing the problem in GDAL master, where
the following query doesn't seem to apply the specified filter:

*ogrinfo "mydata.gdb" -dialect OGRSQL -sql "select * from tablename where
FID=..."*

but this one (without dialect) returns the correct result:

* ogrinfo "mydata.gdb" -sql "select * from tablename where FID=..."  *

I'm a bit uncertain about whether it should be fixed in the generic SQL
layer or the FGDB layer (which doesn't expose OBJECTID in the layer defn.),
since the FGDB layer reports the name "OBJECTID" in GetFIDColumn, therefore
the gen SQL layer could automatically consider it in the layer definition
when processing the attribute filter. When adding the OBJECTID to the query
it fixes the problem, but that seems to be a bit awkward:

*ogrinfo "mydata.gdb" -dialect OGRSQL -sql "select *, OBJECTID from
tablename where FID=..." *

Any ideas in this topic would be helpful?

Best regards,

Tamas


Tamas Szekeres  ezt írta (időpont: 2019. szept. 5.,
Cs, 22:51):

> Hi,
>
> I'm not sure if bForwardWhereToSourceLayer should not be set in this case,
> since the special field FID in pszWHEREIn has already been replaced.
> Or the OpenFileGDB driver should indeed expose OBJECTID as a column
> according to #4253
>
> Best regards,
>
> Tamas
>
>
>  Even Rouault  ezt írta (időpont: 2019.
> szept. 5., Cs, 13:13):
>
>> On mercredi 4 septembre 2019 22:13:40 CEST Tamas Szekeres wrote:
>> > Hi,
>> >
>> > It looks like  the sql queries with -dialect "OGRSQL" doesn't seem to
>> work
>> > as expected. When I specify the FID in the where clause, it doesn't
>> filter
>> > anything. The same query is also described as a solution in the
>> following
>> > ticket https://trac.osgeo.org/gdal/ticket/4253 but I doubt if that
>> works at
>> > all.
>> >
>> > The code causing this problem is fairly generic (ogr_gensql.cpp)
>> >
>> > if( psSelectInfo->where_expr && pszDialect != nullptr &&
>> > EQUAL(pszDialect, "OGRSQL") )
>> > {
>> > int nMinIndexForSpecialField =
>> > poSrcLayer->GetLayerDefn()->GetFieldCount();
>> > bForwardWhereToSourceLayer =
>> > !OGRGenSQLResultsLayerHasSpecialField
>> > (psSelectInfo->where_expr,
>> > nMinIndexForSpecialField);
>> > }
>> > if (bForwardWhereToSourceLayer)
>> > pszWHERE = CPLStrdup(pszWHEREIn);
>> > else
>> > pszWHERE = nullptr;
>> >
>> > In the "where" expression, the FID field is thanslated to OBJECTID and
>> it
>> > is now treated as a special field, therefore the "where" expression is
>> not
>> > being passed to the driver.
>> >
>> > I'm also unsure if that is a correct action to omit passing "where" to
>> the
>> > layer instead of providing an error message.
>>
>> Actually the where isn't completely discarded. It is set on the
>> OGRGenSQLResultsLayer per
>>
>> if( !bForwardWhereToSourceLayer )
>> OGRGenSQLResultsLayer::SetAttributeFilter( pszWHEREIn );
>>
>> around line 492
>>
>> The issue is that the GenSQL layer has no FID column set, and thus this
>> filter
>> fails. One could potentially set the FID Column name on it from the
>> source
>> layer, but that wouldn't be really appropriate in the case of JOIN. That
>> said
>> I see a poDstFeat->SetFID( poSrcFeat->GetFID() ); at line 1332 of
>> TranslateFeature(), so...
>>
>> (there might have been other fixes since #4253 that have made this case
>> broken)
>>
>> --
>> Spatialys - Geospatial professional services
>> http://www.spatialys.com
>>
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] New Warnings from GTIFF output

2020-01-23 Thread Even Rouault 
On jeudi 23 janvier 2020 14:10:03 CET Brian wrote:
> Strange because the files were created by GDAL I created the COG using
> these steps ( original file is too large it is an georaster export). Do you
> see anything that would introduce the corrupt data? What can be done to fix
> it? I attached each file

(Re-adding the list for concluding this thread)

Thanks for the reproducer. The warning was indeed not legit in that situation. 
Just fixed in internal libtiff of GDAL master & libtiff gitlab master

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] New Warnings from GTIFF output

2020-01-23 Thread Even Rouault 
On jeudi 23 janvier 2020 13:26:01 CET Brian wrote:
> Currently using the master and received these warnings when running
> gdal_translate
> 
> "Warning 1: TIFFReadDirectory:Invalid data type for tag TileByteCounts
> Warning 1: TIFFReadDirectory:Invalid data type for tag TileOffsets"
> 
> The exact command was
> gdal_translate -stats raster_trans.tif raster_cog.tif -co TILED=YES -co
> BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co NUM_THREADS=ALL_CPUS -co
> BIGTIFF=IF_SAFER -co COMPRESS=DEFLATE -co COPY_SRC_OVERVIEWS=YES --config
> GDAL_TIFF_OVR_BLOCKSIZE 256
> 
> Is there anyway to suppress warnings and only return errors?

Not that I can think of, at least not without writing a custom error handler.

This warning is a consequence of

commit 5acc54b38bf9bcd45ae8c475f19eb4b8d44ce0bf
Author: Even Rouault 
Date:   Sun Jan 12 14:54:40 2020 +0100

Internal libtiff: _TIFFPartialReadStripArray: add support for non-
conformant SLONG8 data type (fixes #2165)

Such files are non conformant and could be rejected by other software (and for 
GDAL versions >= 3.0 using internal libtiff or libtiff 4.1.0 and before that 
fix, that file couldn't be read at all), so this warning is legit and should 
be taken seriously.

> These warnings did not show up in 2.4.0 and 3.1.0 (is that equal to 3.0.1)?

3.1.0 is not yet released and is the master version, different of 3.0.1

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] New Warnings from GTIFF output

2020-01-23 Thread Brian 
Currently using the master and received these warnings when running
gdal_translate

"Warning 1: TIFFReadDirectory:Invalid data type for tag TileByteCounts
Warning 1: TIFFReadDirectory:Invalid data type for tag TileOffsets"

The exact command was
gdal_translate -stats raster_trans.tif raster_cog.tif -co TILED=YES -co
BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co NUM_THREADS=ALL_CPUS -co
BIGTIFF=IF_SAFER -co COMPRESS=DEFLATE -co COPY_SRC_OVERVIEWS=YES --config
GDAL_TIFF_OVR_BLOCKSIZE 256

Is there anyway to suppress warnings and only return errors?
These warnings did not show up in 2.4.0 and 3.1.0 (is that equal to 3.0.1)?
I can try to send the raw data used if needed.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGR ExecuteSQL

2020-01-23 Thread Mateusz Loskot 
The new docs are searchable

https://gdal.org/search.html?q=ExecuteSQL


Mateusz Loskot, mate...@loskot.net
(Sent from mobile, may suffer from top-posting)

On Thu, 23 Jan 2020, 18:04 Andrew Bell,  wrote:

>
> Thanks.  Just located it in the documentation.  Google was not my friend
> in this case.
>
> Sorry to bother,
>
> On Thu, Jan 23, 2020 at 1:03 PM Even Rouault 
> wrote:
>
>> On jeudi 23 janvier 2020 12:49:38 CET Andrew Bell wrote:
>> > Hi,
>> >
>> > This function returns a pointer to a layer.  Does the dataset retain
>> > ownership of the layer, or does the user need to make sure that the
>> layer
>> > gets deleted?
>> >
>>
>>
>> https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset10ExecuteSQLEPKcP11OGRGeometryPKc
>>
>> """
>> Return
>> an OGRLayer containing the results of the query. Deallocate with
>> ReleaseResultSet().
>> """
>>
>> --
>> Spatialys - Geospatial professional services
>> http://www.spatialys.com
>>
>
>
> --
> Andrew Bell
> andrew.bell...@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] OGR ExecuteSQL

2020-01-23 Thread Andrew Bell 
Thanks.  Just located it in the documentation.  Google was not my friend in
this case.

Sorry to bother,

On Thu, Jan 23, 2020 at 1:03 PM Even Rouault 
wrote:

> On jeudi 23 janvier 2020 12:49:38 CET Andrew Bell wrote:
> > Hi,
> >
> > This function returns a pointer to a layer.  Does the dataset retain
> > ownership of the layer, or does the user need to make sure that the layer
> > gets deleted?
> >
>
>
> https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset10ExecuteSQLEPKcP11OGRGeometryPKc
>
> """
> Return
> an OGRLayer containing the results of the query. Deallocate with
> ReleaseResultSet().
> """
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
>


-- 
Andrew Bell
andrew.bell...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGR ExecuteSQL

2020-01-23 Thread Even Rouault 
On jeudi 23 janvier 2020 12:49:38 CET Andrew Bell wrote:
> Hi,
> 
> This function returns a pointer to a layer.  Does the dataset retain
> ownership of the layer, or does the user need to make sure that the layer
> gets deleted?
> 

https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset10ExecuteSQLEPKcP11OGRGeometryPKc

"""
Return
an OGRLayer containing the results of the query. Deallocate with 
ReleaseResultSet().
"""

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] OGR ExecuteSQL

2020-01-23 Thread Andrew Bell 
Hi,

This function returns a pointer to a layer.  Does the dataset retain
ownership of the layer, or does the user need to make sure that the layer
gets deleted?

Thanks,

-- 
Andrew Bell
andrew.bell...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

2020-01-23 Thread Even Rouault 
On jeudi 23 janvier 2020 18:35:40 CET Edzer Pebesma wrote:
> On 1/23/20 5:58 PM, Even Rouault wrote:
> >> Thanks! Yes, with SRS->importFromUserInput("EPSG:3857"); etc this works
> >> fine; I'm OK with the warning, I just don't see how the subsequent error
> >> relates to syntax that is deprecated.
> > 
> > I'm looking at this, but haven't yet an explanation. This is very subtle,
> > and seems to be linked specifically to EPSG:3857

OK, memory helping, I finally identified the root cause. This is a bug in 
PROJ. Fix submitted per https://github.com/OSGeo/PROJ/pull/1873

While being deterministic, it indeed looks like random in practice. I'm not 
even sure the workaround I proposed would work reliably in all situations with 
calls to GDAL/PROJ that would have been attempted before (due to caching of 
objects inside PROJ context). So either wait for a version with
https://github.com/OSGeo/PROJ/pull/1873 applied, or migrate away from use of 
the deprecated +init=epsg: syntax

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

2020-01-23 Thread Edzer Pebesma 


On 1/23/20 5:58 PM, Even Rouault wrote:
>> Thanks! Yes, with SRS->importFromUserInput("EPSG:3857"); etc this works
>> fine; I'm OK with the warning, I just don't see how the subsequent error
>> relates to syntax that is deprecated.
> 
> I'm looking at this, but haven't yet an explanation. This is very subtle,
> and seems to be linked specifically to EPSG:3857

I don't think so: I see the same problem if I use +init=epsg:28992
rather than +init=epsg:3875.

Your advice of improving this C++ program is well meant, but fixing this
C++ program is not my problem, the problem is how we use GDAL's SRS
interface in R packages. I saw the error occuring at pretty random
occasions, but always errors _after_ one or more +init=epsg: strings
had been used to initialize an SRS.

> 
> I've modified your reproducer as the following:
> 
> int main() {
> 
>OGRSpatialReference *aSRS = new OGRSpatialReference;
> 
>aSRS->importFromWkt("GEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 
> 1984\",ELLIPSOID[\"WGS 
> 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122");
> 
>OGRSpatialReference *bSRS = new OGRSpatialReference;
> #ifdef WORKING
>OGRSpatialReference srs;
>srs.importFromProj4("+init=epsg:3857");
>char* wktb = NULL;
>srs.exportToPrettyWkt();
> #else
>bSRS->importFromProj4("+init=epsg:3857");
>const char* wktb = "PROJCRS[\"WGS 84 / Pseudo-Mercator\",BASEGEOGCRS[\"WGS 
> 84\",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[\"unnamed\",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[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001";
> #endif
>bSRS->importFromWkt((const char *) wktb);
> 
>OGRCoordinateTransformation *ct =
> OGRCreateCoordinateTransformation(aSRS, bSRS);
>if (ct == NULL) {
> printf("ct NULL\n");
> exit(1);
>}
>exit(0);
> }
> 
> So the error is linked to having importFromProj4() and then importFromWkt()
> This doesn't make sense as the later should cancel the former, so there's some
> side effect of importFromProj4() that has later consequences.
> If you define WORKING, which will do the importFromProj4() + 
> exportToPrettyWkt() 
> in a temporary object, and import the resulting WKT in the final bSRS, then
> it works. Doesn't make sense either but could be a workaround
> 
> But Sean's advice is definitely the way forward . Use importFromEPSG() or
> SetFromUserInput("EPSG:")
> And SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) if you don't
> want EPSG compliant axis definitions.
> 
> Even
> 

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081


pEpkey.asc
Description: application/pgp-keys
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

2020-01-23 Thread Even Rouault 
> Thanks! Yes, with SRS->importFromUserInput("EPSG:3857"); etc this works
> fine; I'm OK with the warning, I just don't see how the subsequent error
> relates to syntax that is deprecated.

I'm looking at this, but haven't yet an explanation. This is very subtle,
and seems to be linked specifically to EPSG:3857

I've modified your reproducer as the following:

int main() {

   OGRSpatialReference *aSRS = new OGRSpatialReference;

   aSRS->importFromWkt("GEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 
1984\",ELLIPSOID[\"WGS 
84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122");

   OGRSpatialReference *bSRS = new OGRSpatialReference;
#ifdef WORKING
   OGRSpatialReference srs;
   srs.importFromProj4("+init=epsg:3857");
   char* wktb = NULL;
   srs.exportToPrettyWkt();
#else
   bSRS->importFromProj4("+init=epsg:3857");
   const char* wktb = "PROJCRS[\"WGS 84 / Pseudo-Mercator\",BASEGEOGCRS[\"WGS 
84\",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[\"unnamed\",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[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001";
#endif
   bSRS->importFromWkt((const char *) wktb);

   OGRCoordinateTransformation *ct =
OGRCreateCoordinateTransformation(aSRS, bSRS);
   if (ct == NULL) {
printf("ct NULL\n");
exit(1);
   }
   exit(0);
}

So the error is linked to having importFromProj4() and then importFromWkt()
This doesn't make sense as the later should cancel the former, so there's some
side effect of importFromProj4() that has later consequences.
If you define WORKING, which will do the importFromProj4() + 
exportToPrettyWkt() 
in a temporary object, and import the resulting WKT in the final bSRS, then
it works. Doesn't make sense either but could be a workaround

But Sean's advice is definitely the way forward . Use importFromEPSG() or
SetFromUserInput("EPSG:")
And SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) if you don't
want EPSG compliant axis definitions.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

2020-01-23 Thread Edzer Pebesma 


On 1/23/20 5:29 PM, Sean Gillies wrote:
> Hi,
> 
> On Thu, Jan 23, 2020 at 8:01 AM Edzer Pebesma
> mailto:edzer.pebe...@uni-muenster.de>>
> wrote:
> 
> The combination of GDAL 3.0.3 and PROJ 6.3.0 arriving on debian
> platforms is causing some havoc for some spatial R packages. I could
> bring one particular problem back to the following C++ program:
> 
> #include 
> #include "ogrsf_frmts.h"
> #include 
> 
> int main() {
> 
>    OGRSpatialReference *aSRS = new OGRSpatialReference;
>    OGRSpatialReference *bSRS = new OGRSpatialReference;
> 
>    char *pszSRS = NULL, *wkta = NULL, *wktb = NULL;
> 
>    aSRS->importFromProj4("+init=epsg:4326");
>    aSRS->exportToPrettyWkt();
>    aSRS->importFromWkt((const char *) wkta);
> 
>    bSRS->importFromProj4("+init=epsg:3857");
>    bSRS->exportToPrettyWkt();
>    bSRS->importFromWkt((const char *) wktb);
> 
>    OGRCoordinateTransformation *ct =
> OGRCreateCoordinateTransformation(aSRS, bSRS);
>    if (ct == NULL) {
>         printf("ct NULL\n");
>         exit(1);
>    }
>    exit(0);
> }
> 
> which outputs:
> 
> Warning 1: +init=epsg: syntax is deprecated. It might return a CRS
> with a non-EPSG compliant axis order.
> ERROR 1: PROJ: proj_create_operations: At least one of the operation
> lacks a source and/or target CRS
> ERROR 6: Cannot find coordinate operations from `GEOGCRS["WGS
> 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS
> 
> 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122'
> to `PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",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["unnamed",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001'
> ct NULL
> 
> What am I doing wrong?
> 
> 
> If you pass "EPSG:4326" to importFromProj4 instead of "+init=epsg:4326"
> (and likewise for EPSG:3857) the deprecated syntax warning won't be
> printed. I expect that the error would be cleared up too, at least that
> is my experience with the latest in these libraries.

Thanks! Yes, with SRS->importFromUserInput("EPSG:3857"); etc this works
fine; I'm OK with the warning, I just don't see how the subsequent error
relates to syntax that is deprecated.

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081


pEpkey.asc
Description: application/pgp-keys
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] OPeNDAP support with libdap

2020-01-23 Thread Wilco K 
LS,

which GDAL version includes the DODS driver (libdap)?
We need to query the NOAA GrADS Data Server for GFS data with c#.

Kind regards

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

Re: [gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

2020-01-23 Thread Sean Gillies 
Hi,

On Thu, Jan 23, 2020 at 8:01 AM Edzer Pebesma 
wrote:

> The combination of GDAL 3.0.3 and PROJ 6.3.0 arriving on debian
> platforms is causing some havoc for some spatial R packages. I could
> bring one particular problem back to the following C++ program:
>
> #include 
> #include "ogrsf_frmts.h"
> #include 
>
> int main() {
>
>OGRSpatialReference *aSRS = new OGRSpatialReference;
>OGRSpatialReference *bSRS = new OGRSpatialReference;
>
>char *pszSRS = NULL, *wkta = NULL, *wktb = NULL;
>
>aSRS->importFromProj4("+init=epsg:4326");
>aSRS->exportToPrettyWkt();
>aSRS->importFromWkt((const char *) wkta);
>
>bSRS->importFromProj4("+init=epsg:3857");
>bSRS->exportToPrettyWkt();
>bSRS->importFromWkt((const char *) wktb);
>
>OGRCoordinateTransformation *ct =
> OGRCreateCoordinateTransformation(aSRS, bSRS);
>if (ct == NULL) {
> printf("ct NULL\n");
> exit(1);
>}
>exit(0);
> }
>
> which outputs:
>
> Warning 1: +init=epsg: syntax is deprecated. It might return a CRS
> with a non-EPSG compliant axis order.
> ERROR 1: PROJ: proj_create_operations: At least one of the operation
> lacks a source and/or target CRS
> ERROR 6: Cannot find coordinate operations from `GEOGCRS["WGS
> 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS
>
> 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122'
> to `PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",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["unnamed",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001'
> ct NULL
>
> What am I doing wrong?


If you pass "EPSG:4326" to importFromProj4 instead of "+init=epsg:4326"
(and likewise for EPSG:3857) the deprecated syntax warning won't be
printed. I expect that the error would be cleared up too, at least that is
my experience with the latest in these libraries.

-- 
Sean Gillies
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

2020-01-23 Thread Edzer Pebesma 
The combination of GDAL 3.0.3 and PROJ 6.3.0 arriving on debian
platforms is causing some havoc for some spatial R packages. I could
bring one particular problem back to the following C++ program:

#include 
#include "ogrsf_frmts.h"
#include 

int main() {

   OGRSpatialReference *aSRS = new OGRSpatialReference;
   OGRSpatialReference *bSRS = new OGRSpatialReference;

   char *pszSRS = NULL, *wkta = NULL, *wktb = NULL;

   aSRS->importFromProj4("+init=epsg:4326");
   aSRS->exportToPrettyWkt();
   aSRS->importFromWkt((const char *) wkta);

   bSRS->importFromProj4("+init=epsg:3857");
   bSRS->exportToPrettyWkt();
   bSRS->importFromWkt((const char *) wktb);

   OGRCoordinateTransformation *ct =
OGRCreateCoordinateTransformation(aSRS, bSRS);
   if (ct == NULL) {
printf("ct NULL\n");
exit(1);
   }
   exit(0);
}

which outputs:

Warning 1: +init=epsg: syntax is deprecated. It might return a CRS
with a non-EPSG compliant axis order.
ERROR 1: PROJ: proj_create_operations: At least one of the operation
lacks a source and/or target CRS
ERROR 6: Cannot find coordinate operations from `GEOGCRS["WGS
84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS
84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122'
to `PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",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["unnamed",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["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001'
ct NULL

What am I doing wrong?
-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081


pEpkey.asc
Description: application/pgp-keys
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Querying many raster pixels efficiently from Python

2020-01-23 Thread Daniel Evans 
Hi,

I have a large (global, 30m resolution, 50GB+) GeoTIFF dataset, from which I 
need to read many (millions) of pixel values at given input coordinates. I've 
got reasonable performance out of the code, about a million queries over five 
minutes, but:


  1.  There are actually twelve separate datasets of this size to query, not 
just one, so it takes approximately an hour.
  2.  This is by far the slowest portion of the program, and the users demand 
speed!
  3.  The users would also like to move towards higher resolution datasets, 
which we see run about 5x slower.
  4.  When querying the data on a particular piece of network storage mounted 
as part of the local filesystem, we see a slowdown approaching two orders of 
magnitude - bulk file copies off the network storage are reasonable, but each 
IO request shows a significant overhead (up to a second), and GDAL is sending 
one for each coordinate queried.

The implementation is in Python, directly calling down to GDAL. The short, 
long-running snippet of code which performs the actual queries the dataset, 
having converted real-world coordinates to pixels, is:

value_arrays = (
raster_ds.ReadAsArray(
xoff=coord[0] - buffer_size,
yoff=coord[1] - buffer_size,
xsize=npix,
ysize=npix
) for coord in offsets
)

There are a few things that are probably worth noting:

  1.  It is not necessarily a single pixel that is being read - for each 
coordinate, the program may be asked to get all pixel values within a given 
radius (typically a couple of pixels), and use some function to summarise these 
into a single value (mean, median, ...). GDAL currently returns a numpy array 
for each query, which is passed to the user-specified function after the 
snippet above.
  2.  The dataset is made up of 2048x2048 LZW-Compressed tiles containing 
floating point data (essentially conforming to COG, but with no overviews), 
grouped together in a VRT (performance is identical with plain GeoTIFFs, 
though).

Multiprocessing has not been found to help - we actually lose throughput as the 
disk read head is moving back and forth constantly. Better hardware (especially 
SSDs) is known to help, but no one wants to pay for that.

We see no particular performance difference from setting 
GDAL_DISABLE_READDIR_ON_OPEN=TRUE, and GDAL_CACHEMAX is left at the default 5% 
(64GB+ RAM available).



Does the Python interface to GDAL provide a way to supply a large number of 
offsets and get blocks of pixels back, avoiding the need to come back up to 
Python after each query? (I suspect not)
Is there some way to optimise GDAL so that queries of files on the mounted 
network storage are more efficient?



Dr. Daniel Evans
Software Developer

Skype


T +44 (0) 1756 799919
www.jbarisk.com

[Visit our website] 
[http://www.jbagroup.co.uk/imgstore/JBA-Email-Sig-Icons-LINKEDIN.png]  [Follow 
us on Twitter] 
Our postal address and registered office is JBA Risk Management Limited, 1 
Broughton Park, Old Lane North, Broughton, Skipton, North Yorkshire BD23 3FD.

Find out more about us here: www.jbarisk.com and 
follow us on Twitter @JBARisk and 
LinkedIn

The JBA Group supports the JBA Trust.

All JBA Risk Management's email messages contain confidential information and 
are intended only for the individual(s) named. If you are not the named 
addressee you should not disseminate, distribute or copy this e-mail.
Please notify the sender immediately by email if you have received this email 
by mistake and delete this email from your system.


JBA Risk Management Limited is registered in England, company number 07732946, 
1 Broughton Park, Old Lane North, Broughton, Skipton, North Yorkshire, BD23 
3FD, Telephone: +441756799919


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

Re: [gdal-dev] "gdaladdo: tif_dirwrite.c:2449: TIFFWriteDirectoryTagData: Assertion `datalength<0x80000000UL' failed."

2020-01-23 Thread jratike80 
Hi,

Do I estimate right that you have about 125 terabytes image data as
uncompressed , and the first level overview would make about 32 terabytes?
As JPEG compressed it would still make about 3 terabyte .ovr file. That is
not too much for BigTIFF that has a size limit at 18000 petabytes, but it
may still be unpractical.

I would rather split the source data into a bit smaller VRT mosaics with
size around 100x100 km and start creating overviews for those. Or create a
few first levels of overviews for the source tiffs and then materialize
100x100 km areas into geotiffs havin pixel size of 16 meters with
gdal_translate. Or something else, but 3 TB overview file would not be my
first goal. Of course it can be the best alternative for your exact use
case.

-Jukka Rahkonen-


mapmonster wrote
> Hi,
> 
> I'm trying to build overviews on a large vrt file, but I'm getting the
> following error message
> "gdaladdo: tif_dirwrite.c:2449: TIFFWriteDirectoryTagData: Assertion
> `datalength<0x8000UL' failed."
> 
> Size is 3812500, 8203125





--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] "gdaladdo: tif_dirwrite.c:2449: TIFFWriteDirectoryTagData: Assertion `datalength<0x80000000UL' failed."

2020-01-23 Thread Even Rouault 
Mats,

> I'm trying to build overviews on a large vrt file, but I'm getting the > 
following error message
>"gdaladdo: tif_dirwrite.c:2449: TIFFWriteDirectoryTagData: Assertion > 
`datalength<0x8000UL' failed."

Ouch, you hit a bug in libtiff. What you are trying to create is an absolutely 
large BigTIFF file.

Given 128x128 tiles and your raster dimensions, the size of each of the 
TileByteCount/TileOffset TIFF array would be

(3812500 / 128.) * (8203125 / 128.) * sizeof(double) = 15.2 GB !

For tag data, libtiff has a limitation to 2 GB (this is a - reasonable - 
implementation limitation, the format itself could allow to 2^63). Ideally, it 
should have rejected from the start this attempt instead of crashing. You may 
want to file a ticket about that.

You can add --config GDAL_TIFF_OVR_BLOCKSIZE 1024 (or a larger value) to 
increase the size of overview tiles, and thus reduce the size of each array to 
a more reasonable ~ 230 MB.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] "gdaladdo: tif_dirwrite.c:2449: TIFFWriteDirectoryTagData: Assertion `datalength<0x80000000UL' failed."

2020-01-23 Thread Mats Högström 
Hi,

I'm trying to build overviews on a large vrt file, but I'm getting the 
following error message
"gdaladdo: tif_dirwrite.c:2449: TIFFWriteDirectoryTagData: Assertion 
`datalength<0x8000UL' failed."

Operating system ubuntu server 18.04
256 Gb memory
GDAL version 2.4

This is what I do:
- create a text file with all the file names
- gdalbuildvrt mosaic.vrt --optfile 
- gdaladdo -q -ro --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW 
RGB --config INTERLEAVE_OVERVIEW PIXEL --config BIGTIFF_OVERVIEW YES mosaic.vrt 
2

The vrt file references about 15000 4-band tiffs with 16 cm resolution.
gdalinfo mosaic.vrt
Driver: VRT/Virtual Raster
Files: 
/mnt/data/lm/orto4band016/raster/epsg3006/2019/delivery/ortofoto/mosaic.vrt
   
/mnt/data/lm/orto4band016/raster/epsg3006/2019/delivery/ortofoto/70_4/700_47_0075_2019.tif
|
|
|
   
/mnt/data/lm/orto4band016/raster/epsg3006/2019/delivery/ortofoto/72_8/728_80_0050_2019.tif
Size is 3812500, 8203125
Coordinate System is:
PROJCS["SWEREF99 TM",
GEOGCS["SWEREF99",
DATUM["SWEREF99",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6619"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4619"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",50],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3006"]]
Origin = (312500.000,754.000)
Pixel Size = (0.166,-0.166)
Corner Coordinates:
Upper Left  (  312500.000, 754.000) ( 10d31'41.87"E, 67d54'48.45"N)
Lower Left  (  312500.000, 6227500.000) ( 11d58'51.51"E, 56d 9'20.09"N)
Upper Right (  922500.000, 754.000) ( 24d59'59.48"E, 67d40' 2.67"N)
Lower Right (  922500.000, 6227500.000) ( 21d46'52.71"E, 56d 0'21.07"N)
Center  (  617500.000, 6883750.000) ( 17d14'55.13"E, 62d 4' 3.41"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  NoData Value=0
Band 4 Block=128x128 Type=Byte, ColorInterp=Undefined
  NoData Value=0

Any ideas why this happens?

Best regards
/mats
---
När du skickar e-post till SLU så innebär detta att SLU behandlar dina 
personuppgifter. För att läsa mer om hur detta går till, klicka här 

E-mailing SLU will result in SLU processing your personal data. For more 
information on how this is done, click here 

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

[gdal-dev] gdalinfo not working on Mageia

2020-01-23 Thread David GEIGER 
Hi,

On Mageia we have an issue with gdalinfo (and some others scripts), during
build it creates a bash script instead of a real binary, I checked what is
going wrong but without success:

$ gdalinfo --version
/usr/bin/gdalinfo: error: '/usr/bin/.libs/gdalinfo' does not exist

Mageia7 and Mageia Cauldron
gdal 2.4.3

https://bugs.mageia.org/show_bug.cgi?id=25809#c3

Here our spec file:

http://svnweb.mageia.org/packages/cauldron/gdal/current/SPECS/gdal.spec?revision=1478884=markup



Regards,

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