Re: [gdal-dev] [EXT] Re: Unknown projection

[Matt.Wilkie via gdal-dev]

Thanks Even!


Matt




From: Even Rouault 
Sent: Wednesday, January 31, 2024 12:54
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org 

Subject: [EXT] Re: [gdal-dev] Unknown projection


Matt,

surprisingly or not, the GTiff driver doesn't try to get the SRS from a .prj 
side car file.

So just do

gdal_edit.py 047_0001.tif -a_srs "$(cat 047_0001.prj)"

Even

Le 31/01/2024 à 20:22, Matt.Wilkie via gdal-dev a écrit :

Hello gdal,



I have a collection of aerial photos for which neither gdal nor arc pro 
recognize the projection. To my eye the .prj files look ok. Attached txt is 
gdalinfo v3.8.0 report along with the .prj. I manually formatted the below for 
readability:



PROJCS["UTM_07N",

GEOGCS["NAD83-CSRS Epoch 2006 (7 Param+Vel)",

DATUM["NAD83-CSRS Epoch 2006 (7 Param+Vel)",

SPHEROID["GRS 80", 6378137,298.2572221010002],

TOWGS84[1.0026,-1.9083,-0.5165, 1.288877171229702e-07, 
8.998141921392988e-09, 5.376098909823645e-08, -1.18e-09]],

PRIMEM["Greenwich",0],

UNIT["degree",0.0174532925199433]

],

PROJECTION["Transverse_Mercator"],

PARAMETER["latitude_of_origin",0],

PARAMETER["central_meridian",-141],

PARAMETER["scale_factor",0.9996],

PARAMETER["false_easting",50],

PARAMETER["false_northing",0],

UNIT["metres",1]

]



Thanks!

Matt Wilkie

Geomatics Developer & Administrator

Environment | Technology, Innovation and Mapping

T 867-667-8133 | Yukon.ca<http://yukon.ca/>

Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.





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


--
http://www.spatialys.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

[gdal-dev] Unknown projection

[Matt.Wilkie via gdal-dev]

Hello gdal,

I have a collection of aerial photos for which neither gdal nor arc pro 
recognize the projection. To my eye the .prj files look ok. Attached txt is 
gdalinfo v3.8.0 report along with the .prj. I manually formatted the below for 
readability:

PROJCS["UTM_07N",
GEOGCS["NAD83-CSRS Epoch 2006 (7 Param+Vel)",
DATUM["NAD83-CSRS Epoch 2006 (7 Param+Vel)",
SPHEROID["GRS 80", 6378137,298.2572221010002],
TOWGS84[1.0026,-1.9083,-0.5165, 1.288877171229702e-07, 
8.998141921392988e-09, 5.376098909823645e-08, -1.18e-09]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]
],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-141],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",50],
PARAMETER["false_northing",0],
UNIT["metres",1]
]

Thanks!
Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.



047_0001.prj
Description: 047_0001.prj
Driver: GTiff/GeoTIFF
Files: 047_0001.tif
   047_0001.tfw
Size is 14592, 25728
GeoTransform =
  551461.43038, -0.2658, -0.00426
  6961619.656405, -0.00103, 0.2651
Metadata:
  TIFFTAG_MAXSAMPLEVALUE=255
  TIFFTAG_MINSAMPLEVALUE=0
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  551461.430, 6961619.656) 
Lower Left  (  551351.829, 6968440.149) 
Upper Right (  547582.877, 6961604.627) 
Lower Right (  547473.276, 6968425.119) 
Center  (  549467.353, 6965022.388) 
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [EXT] Happy Birthday GDAL!

[Matt.Wilkie via gdal-dev]

Oh wow. Happy birthday gdal! And thank you Frank Warmerdam.

That means I've been using gdal for just over 20 years, though the oldest 
extant public record I could find is 2008 
(ref). 
This makes me feel old, in a good way.

Matt
Geomatics Developer and Administrator | Environment | T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

From: gdal-dev  On Behalf Of Javier Jimenez 
Shaw via gdal-dev
Sent: Friday, October 20, 2023 1:01 AM
To: gdal dev 
Subject: [EXT] [gdal-dev] Happy Birthday GDAL!

GDAL got 25 years old on 2023-10-17 !
I found it in this post 
https://geoobserver.wordpress.com/2023/10/17/25-jahre-gdal-happy-birthday/
.___ ._ ..._ .. . ._.  .___ .. __ . _. . __..  ...  ._ .__
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] JPEG compressed COGs nodata

[Matt.Wilkie]

For those not wanting or needing to take the processing hit of running 
nearblack, use gdalwarp and -dstalpha to create the alpha channel, and use 
intermediate VRT to save space and more time. Something like this:

gdalwarp -srcnodata 0 -dstalpha -of vrt source.tif xx-interim.vrt

gdal_translate -of COG {...} xx-interim.vrt final.tif


-Matt

From: gdal-dev  On Behalf Of Travis Kirstine
Sent: July 22, 2022 9:04 AM
To: gdal dev 
Subject: Re: [gdal-dev] JPEG compressed COGs nodata


You don't often get email from 
traviskirst...@gmail.com. Learn why this is 
important

To answer my own question I believe the answer is that you need to create an 
alpha band when working with RGB datasets to have JPEG compressed COGs without 
compression artifacts in the nodata areas.  Creating a internal mask on the 
source data may work as well but I haven't tried.

# add a alpha band using the nearblack utility
nearblack -setalpha rgb_nodata.tif -o rgba.tif
# create the COG
gdal_translate -of COG  -co COMPRESS=JPEG -co 
TILING_SCHEME=GoogleMapsCompatible -co NUM_THREADS=ALL_CPUS -co BIGTIFF=YES 
rgba.tif cog.tif

Regards




On Tue, 19 Jul 2022 at 15:19, Travis Kirstine 
mailto:traviskirst...@gmail.com>> wrote:
I've been trying to create a COG with JPEG compression with transparent nodata 
values without much success.  Is this possible without creating a secondary 
mask?

Any hints?

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

Re: [gdal-dev] Linux - bash script - Sentinel 1 processing - dark images

[Matt.Wilkie]

> Does anybody have some hints for me how to improve the exposure/contrast of 
> the images?

This doesn't help directly with constructing any command lines you might use, 
but I found the "Gentle Introduction to GDAL" series by Robert Simmon very 
helpful in determining what to do.

Here's Part 5, Saturation and Sharpness:
https://medium.com/@robsimmon/making-sense-of-satellite-data-part-5-saturation-sharpness-e21fac6272c7

In my memory I have a similarly written article on gdal colour balancing that 
did have some command line examples but I can't relocate it at the moment. It 
annoys me that I can't reference it so I'll keep looking. ;-)

cheers,

-Matt
___
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?

[Matt.Wilkie]

If you decide to add an elevation metadata keyword for things you generate let 
us know and I at least will follow suit. ;-)

-Matt

-Original Message-
From: gdal-dev  On Behalf Of Nyall Dawson
Sent: April 29, 2022 6:11 PM
To: gdal dev 
Subject: [gdal-dev] Hints that a raster band represents a DEM?

[You don't often get email from nyall.daw...@gmail.com. Learn why this is 
important at https://aka.ms/LearnAboutSenderIdentification.]

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

Re: [gdal-dev] Standardize gdal-utils scripts return code for "no arguments"

[Matt.Wilkie]

Thank you for the study path Idan, much appreciated and I will use it.

> Bottom line, returncode 2 is from the very last line in argparse.py,
> the builtin ArgumentParser (function ArgumentParser.error()):
> `self.exit(2, _('%(prog)s: error: %(message)s\n') % args)`

Hmmm. Since return 2 is the default for argparse then it seems to we should do 
the same for all the python utils?

-Matt

From: Idan Miara  On Behalf Of Idan Miara
Sent: April 27, 2022 6:24 AM
To: Matt.Wilkie 
Cc: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Standardize gdal-utils scripts return code for "no 
arguments"


You don't often get email from i...@miara.com<mailto:i...@miara.com>. Learn why 
this is important<https://aka.ms/LearnAboutSenderIdentification>

Hi Matt,

Since you mentioned that you don't understand how classes work (i.e. super 
function)
I would highly recommend spending some time to go over some basic Python OOP, 
it's not complicated and I'm confident that it would be very useful.
From a quick search, maybe this link would be useful 
https://realpython.com/python3-object-oriented-programming/<https://can01.safelinks.protection.outlook.com/?url=https%3A%2F%2Frealpython.com%2Fpython3-object-oriented-programming%2F=05%7C01%7CMatt.Wilkie%40yukon.ca%7C44e0362e122546ade38008da28514d88%7C98f515313973490abb70195aa264a2bc%7C0%7C0%7C637866628546448376%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=lldjRlTrTEnWnk4lhViopt%2B2v4cmk1li2JXWDGySfO0%3D=0>
 but there are many other resources.
super() is mainly used in Python as a way to call the superclass' function when 
that you subclass 
(https://realpython.com/python-super/<https://can01.safelinks.protection.outlook.com/?url=https%3A%2F%2Frealpython.com%2Fpython-super%2F=05%7C01%7CMatt.Wilkie%40yukon.ca%7C44e0362e122546ade38008da28514d88%7C98f515313973490abb70195aa264a2bc%7C0%7C0%7C637866628546448376%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=KZ%2B4PamCCg0v1nA2hFrq1LsSeBXAN7oKCKWgLpyW%2FKU%3D=0>)
For example, in `GDALArgumentParser.__init__()` I call `super().__init__` in 
order to access ArgumentParser's constructor (__init__).

GDALArgumentParser is defined as `class 
GDALArgumentParser(argparse.ArgumentParser)`.
So GDALArgumentParser is a superclass of the built in class 
argparse.ArgumentParser.
The objective of this class is mainly to standardize our usage of 
ArgumentParser by adding some functionality that was implemented in many 
different utils on top of ArgumentParser.
As GDALArgumentParser is a thin superclass of ArgumentParser, If you want to 
better understand GDALArgumentParser I advise you to first read the docs here 
https://docs.python.org/3/library/argparse.html<https://can01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdocs.python.org%2F3%2Flibrary%2Fargparse.html=05%7C01%7CMatt.Wilkie%40yukon.ca%7C44e0362e122546ade38008da28514d88%7C98f515313973490abb70195aa264a2bc%7C0%7C0%7C637866628546448376%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=rbnjs2pgbJh5nIb8iYJE2BJQX0v65yEsopaf3iWaZ5g%3D=0>,
 it helped me very much when I coded GDALArgumentParser.

Bottom line, returncode 2 is from the very last line in argparse.py, the 
builtin ArgumentParser (function ArgumentParser.error()):
`self.exit(2, _('%(prog)s: error: %(message)s\n') % args)`
When ArgumentParser encounters an error (like calling it without arguments when 
arguments were required) it calls this function and exits with errorcode 2.
You could see this in action if you copy-paste and run without parameters the 
first example in the ArgumentParser docs 
(https://docs.python.org/3/library/argparse.html<https://can01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdocs.python.org%2F3%2Flibrary%2Fargparse.html=05%7C01%7CMatt.Wilkie%40yukon.ca%7C44e0362e122546ade38008da28514d88%7C98f515313973490abb70195aa264a2bc%7C0%7C0%7C637866628546448376%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=rbnjs2pgbJh5nIb8iYJE2BJQX0v65yEsopaf3iWaZ5g%3D=0>).

The way I found this was just using the debugger in PyCharm and stepping in 
until I found out what's exits the run
(If an exception is raised it makes it much easier as you can easily follow the 
call stack)

pyi are stub files 
(https://peps.python.org/pep-0484/<https://can01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpeps.python.org%2Fpep-0484%2F=05%7C01%7CMatt.Wilkie%40yukon.ca%7C44e0362e122546ade38008da28514d88%7C98f515313973490abb70195aa264a2bc%7C0%7C0%7C637866628546448376%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=YgIYRRFyd1rKQnLQukduYh71LiPtV5XltFow%2B31%2BdkU%3D=0>)
 and you probably encountered them when digging into Python's standard library 
which is not written in pure Python

Re: [gdal-dev] Standardize gdal-utils scripts return code for "no arguments"

[Matt.Wilkie]

Here is the path I'm following in trying to deconstruct GDALArgumentParser.
With file = pc2rgb.py:
r = subprocess.run([sys.executable,
file],
shell=True,
capture_output=True,
text=True,
)
if debug:
print(f'returncode: {r.returncode}')
I get:
returncode: 2
Digging in to pct2rgb.py I see main() calling PCT2RGB class, calling pct2rgb 
function, which should error out with an exception.:
def main(argv=sys.argv):
return PCT2RGB().main(argv)
def doit(self, **kwargs):
return pct2rgb(**kwargs)
def doit(**kwargs):
try:
ds = pct2rgb(**kwargs)
return ds, 0
except:
return None, 1
def pct2rgb(src_filename: PathLikeOrStr, pct_filename: Optional[PathLikeOrStr], 
dst_filename: PathLikeOrStr,
band_number: int = 1, out_bands: int = 3, driver_name: 
Optional[str] = None):
# Open source file
src_ds = open_ds(src_filename)
if src_ds is None:
raise Exception(f'Unable to open {src_filename} ')
...
So I'm thinking that while PCT2RGB class is invoking GDALArgumentParser it's 
exiting early somehow with this 2 return code.
def get_parser(self, argv) -> GDALArgumentParser:
...
return parser
I look at auxillary/gdal_argparse.py: class 
GDALArgumentParser(argparse.ArgumentParser) and I see... well a bunch of stuff 
I only apprehend the shadowy outlines of. Looking at the return statements 
doesn't add light. (It doesn't help that I don't understand how classes work. 
There's this thing called super() but it's not defined anywhere. Shouldn't that 
cause an error? Oh it's from stdlib builtins.pyi. What the heck is a pyi? 
"...It is pitch black. You are likely to be eaten by a grue.")
> findstr "return" ..\auxiliary\gdal_argparse.py
return super().parse_args(args=args, **kwargs)
return shlex.split(arg_line, comments=True)
return self._parser
return kwargs
return kwargs
return 0
return 1
return epilog or None
So my question is: where is returncode: 2 coming from?

-Matt

From: Matt.Wilkie
Sent: April 26, 2022 1:03 PM
To: 'gdal-dev@lists.osgeo.org' 
Subject: RE: Standardize gdal-utils scripts return code for "no arguments"

Hi Folks,

I've converted all the scripts that were using -1 to 1. However when I started 
looking at the ones returning 2 it became less clear what to do. Excepting 
gdal2tiles all of them are using GDALArgumentParser and I don't see where the 
return value is being set. Thinking this might mean the mainline utils might be 
using 2 for no args I checked gdal_translate and gdalwarp, but no, both of 
those use 1.

The scripts that return 2 are:

osgeo_utils\gdal2xyz.py
osgeo_utils\gdal_calc.py
osgeo_utils\gdal_fillnodata.py
osgeo_utils\gdal_polygonize.py
osgeo_utils\pct2rgb.py
osgeo_utils\rgb2pct.py
osgeo_utils\samples\gdallocationinfo.py


I haven't been able to figure out the GDALArgumentParser  parser class so I 
embarked on making everything else use 2 instead of 1. I'm questioning the 
wisdom of that at the moment since so many files are touched. However the 
process has forced me to look more closely at the many scripts and start to 
internalize the various patterns they use. This has been worthwhile even if the 
approach might get abandoned. The "make everything return 2" effort is in 
branch patch-5561-ret2<https://github.com/maphew/gdal/tree/patch-5561-ret2>.

The next message will have more details on how I've tried and failed to 
understand GDALArgumentParser.

-Matt

From: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> On 
Behalf Of matt.wil...@yukon.ca<mailto:matt.wil...@yukon.ca>
Sent: April 4, 2022 4:07 PM
To: gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
Subject: [gdal-dev] Standardize gdal-utils scripts return code for "no 
arguments"

Hi folks,

I'm working on "Standardize gdal-utils scripts return codes #5561" for all the 
scripts in swig/python/gdal-utils. Currently the scripts do not return the same 
status code for "was run without arguments".

It would be good for the same code meant the same thing across all the scripts 
in the package. Given that most of the scripts use `1` now, and that this is in 
line with sys,exit() docs I think it makes sense to make 1 the new standard.

At present we have (return_code, num scripts with that code):

0: 30
1: 56
2:  8
-1:  9

-1 is a special case. In the script code it's written as `return -1` but 
subprocess.run() captures it as `4294967295`. Another special case is 
`gdal_auth.py` sample which with no arguments spawns a web authentication page 
in browser.

Changing the return code will mean anyone who is relying on those in their own 
scripts or programs will need to adjust. This seems to be a small price for the 
gain in harmonization across the utilities, to me. Your thoughts?


[0]: 
https

Re: [gdal-dev] Standardize gdal-utils scripts return code for "no arguments"

[Matt.Wilkie]

Hi Folks,

I've converted all the scripts that were using -1 to 1. However when I started 
looking at the ones returning 2 it became less clear what to do. Excepting 
gdal2tiles all of them are using GDALArgumentParser and I don't see where the 
return value is being set. Thinking this might mean the mainline utils might be 
using 2 for no args I checked gdal_translate and gdalwarp, but no, both of 
those use 1.

The scripts that return 2 are:

osgeo_utils\gdal2xyz.py
osgeo_utils\gdal_calc.py
osgeo_utils\gdal_fillnodata.py
osgeo_utils\gdal_polygonize.py
osgeo_utils\pct2rgb.py
osgeo_utils\rgb2pct.py
osgeo_utils\samples\gdallocationinfo.py


I haven't been able to figure out the GDALArgumentParser  parser class so I 
embarked on making everything else use 2 instead of 1. I'm questioning the 
wisdom of that at the moment since so many files are touched. However the 
process has forced me to look more closely at the many scripts and start to 
internalize the various patterns they use. This has been worthwhile even if the 
approach might get abandoned. The "make everything return 2" effort is in 
branch patch-5561-ret2.

The next message will have more details on how I've tried and failed to 
understand GDALArgumentParser.

-Matt

From: gdal-dev  On Behalf Of 
matt.wil...@yukon.ca
Sent: April 4, 2022 4:07 PM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Standardize gdal-utils scripts return code for "no 
arguments"

Hi folks,

I'm working on "Standardize gdal-utils scripts return codes #5561" for all the 
scripts in swig/python/gdal-utils. Currently the scripts do not return the same 
status code for "was run without arguments".

It would be good for the same code meant the same thing across all the scripts 
in the package. Given that most of the scripts use `1` now, and that this is in 
line with sys,exit() docs I think it makes sense to make 1 the new standard.

At present we have (return_code, num scripts with that code):

0: 30
1: 56
2:  8
-1:  9

-1 is a special case. In the script code it's written as `return -1` but 
subprocess.run() captures it as `4294967295`. Another special case is 
`gdal_auth.py` sample which with no arguments spawns a web authentication page 
in browser.

Changing the return code will mean anyone who is relying on those in their own 
scripts or programs will need to adjust. This seems to be a small price for the 
gain in harmonization across the utilities, to me. Your thoughts?


[0]: 
https://github.com/OSGeo/gdal/issues/5561
[1]: 
https://docs.python.org/3/library/sys.html#sys.exit,
 "Unix programs generally use 2 for command line syntax errors and 1 for all 
other kind of errors"

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | 
Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Convert to min containing bit depth?

[Matt.Wilkie]

Thank you both for this, very helpful.

-Matt

From: Idan Miara 
Sent: April 21, 2022 11:18 PM
To: Mike Taves 
Cc: Matt.Wilkie ; gdal dev 
Subject: Re: [gdal-dev] Convert to min containing bit depth?

Computing the min value is also requited if you have negative values and could 
also be useful If you wanted to optimize that method further by utilising the 
offset parameter (and or the scale, but computing the right combination for an 
optimize lossless compression could be more expensive).

On Fri, 22 Apr 2022, 06:05 Mike Taves, 
mailto:mwto...@gmail.com>> wrote:
On Fri, 22 Apr 2022 at 07:05, 
mailto:matt.wil...@yukon.ca>> wrote:
>
> Idea for a small but useful python tool: scan image for min/max values and 
> convert to smallest possible bit depth without losing values. Surely someone 
> has done something like this already. Any suggestions for where to look for 
> prior art?

This is driver-specific, as certain formats expect multiples of 2
(e.g.) NBITS=1/2/4. But for GTiff, what I typically use in a script is
to find the maximum value, then use "ceil(log(maxval, 2))" to get the
number of bits, e.g.:

from math import ceil, log
from osgeo import gdal

maxval = 17  # for example
nbits = ceil(log(maxval, 2))  # 5

drv = gdal.GetDriverByName("GTiff")
opts = [f"NBITS={maxval}"]
ds = drv.Create(fname, nx, ny, 1, gdal.GDT_Byte, opts)
...

similar can be done with rasterio, passing the keyword
"rasterio.open(fname, 'w', ..., nbits=nbits)"

For the drivers that expect NBITS as a multiple of 2:

nbits = 2**ceil(log(nbits, 2))

If nbits is greater than 8, then UInt16 or UInt32 may be required, as
supported by the driver.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
https://lists.osgeo.org/mailman/listinfo/gdal-d<https://can01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev=05%7C01%7CMatt.Wilkie%40yukon.ca%7C530d4395252944fbff3a08da2427e83d%7C98f515313973490abb70195aa264a2bc%7C0%7C0%7C637862051196803406%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=Wxzbgz3TEIvRtrCiUvFSGc5CUoGIbeeCC8VXt1wtfuA%3D=0>e
 vroom
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Convert to min containing bit depth?

[Matt.Wilkie]

Idea for a small but useful python tool: scan image for min/max values and 
convert to smallest possible bit depth without losing values. Surely someone 
has done something like this already. Any suggestions for where to look for 
prior art?

thanks.
Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] Standardize gdal-utils scripts return code for "no arguments"

[Matt.Wilkie]

Hi folks,

I'm working on "Standardize gdal-utils scripts return codes #5561" for all the 
scripts in swig/python/gdal-utils. Currently the scripts do not return the same 
status code for "was run without arguments".

It would be good for the same code meant the same thing across all the scripts 
in the package. Given that most of the scripts use `1` now, and that this is in 
line with sys,exit() docs I think it makes sense to make 1 the new standard.

At present we have (return_code, num scripts with that code):

0: 30
1: 56
2:  8
-1:  9

-1 is a special case. In the script code it's written as `return -1` but 
subprocess.run() captures it as `4294967295`. Another special case is 
`gdal_auth.py` sample which with no arguments spawns a web authentication page 
in browser.

Changing the return code will mean anyone who is relying on those in their own 
scripts or programs will need to adjust. This seems to be a small price for the 
gain in harmonization across the utilities, to me. Your thoughts?


[0]: https://github.com/OSGeo/gdal/issues/5561
[1]: https://docs.python.org/3/library/sys.html#sys.exit, "Unix programs 
generally use 2 for command line syntax errors and 1 for all other kind of 
errors"

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.


gdal-utils-return-codes.csv
Description: gdal-utils-return-codes.csv
import subprocess
# import glob
import sys
from pathlib import Path

here = r"D:\code\public\gdal\swig\python\gdal-utils\osgeo_utils"
scripts = list(Path(here).glob("**/*.py" ))

for s in scripts:
file = Path(s)
if not 'gdal_auth.py' in file.name:
r = subprocess.run([sys.executable,
file],
shell=True,
capture_output=True,
text=True,
)

print(r.returncode, file)


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

[gdal-dev] Gitpod workspace with Ubuntu GDAL

[Matt.Wilkie]

Gitpod workspace with Ubuntu GDAL

How to open a Gitpod machine with the version of GDAL in the 
Ubuntu software library.

  1.  Point browser at https://gitpod.io/#https://github.com/OSGeo/gdal
  2.  Login with Github credentials
  3.  In a couple of minutes the gitpod linux machine will be fired up and 
ready to use with Visual Studio Code open. There might be a cmake error, ignore 
it for now.
  4.  In the Terminal window install gdal (and whatever else you need).

# install gdal binaries, and git

sudo add-apt-repository -y ppa:ubuntugis/ubuntugis-unstable

sudo apt-get update

sudo apt-get install gdal-bin libgdal-dev



# verify gdal working

gdalinfo --version



# install gdal-utilties

python -m pip install swig/python/gdal-utils



# test gdal-utils working

gdal_edit

Optional: run tests to verify what parts of gdal are working

# Prep

python -m pip install -r requirements.txt



# Ensure gdal source code tree and the installed gdal binaries are the same 
version

git checkout v`gdal-config --version`

# and that python bindings match also

python -c "from osgeo import gdal;print(gdal.__version__)"



# Test only python utilities:

pytest pyscripts/



# All tests:

pytest

Pytest not found? Use python -m pytest instead. In next opening of bash shell 
it will work.



Result of investigation in https://github.com/OSGeo/gdal/issues/5442


Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] Guidelines for optimal image dimensions

[Matt.Wilkie]

Does anyone have guidelines of what the general purpose optimal image 
dimensions in pixels are for any given type? (3 band 8bit RGB, 1 band 16bit 
signed, ...)

At some point the time it takes to do anything in memory on a whole image 
exceeds the time it takes to read the next chunk from storage and repeat. I'm 
curious how to figure out what the point of diminishing returns is. The quest 
behind the question is determining a good default tile size for our image 
repository.

There's a hundred variables that's going to make it a different answer for 
every organization, and that it will change over time as hardware and software 
changes. My hunch is there's a close enough approximation that works for most 
situations much of the time.  Before I embark on  series of trials to figure 
out what our local optimum is I figured I'd ask and see what other work in this 
area there might be.

Your thoughts?
Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] GDAL Python style guide?

[Matt.Wilkie]

Hi All,

What's the closest thing to a style guide for GDAL Python code? I found "RFC 
69: C/C++ Code Formatting"[0] but haven't discovered anything for python yet. 
I'm wondering about preferred line length, autopep8 vs black vs ... formatting 
tools, variable naming patterns, etc.

[0] https://gdal.org/development/rfc/rfc69_cplusplus_formatting.html

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] gdal-utils: About GDAL Python Utilities

[Matt.Wilkie]

I revised according to your feedback and created an RST document from it. I’m 
not too experienced with the syntax so there may be some changes needed. Here’s 
the PR:
https://github.com/OSGeo/gdal/pull/5222

-Matt

From: Idan Miara 
Sent: February 1, 2022 12:30 PM
To: Matt.Wilkie 
Cc: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] gdal-utils: About GDAL Python Utilities


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



Hi Matt,

Thanks for looking into this.

Yes, you got it right.
You can find some documentation in the following pages, though it certainly 
could be approved as per your suggestions, maybe also by using some text from 
the RFC.
https://gdal.org/api/index.html#python-api
https://gdal.org/programs/index.html#programs (note the programs that end with 
.py)
https://gdal.org/api/python_samples.html
These doc pages are auto generated from the sources here:
https://github.com/OSGeo/gdal/tree/master/doc/source

> Osgeo_utils doesn’t have a top level document yet.
I'm not sure what you mean by a top level document (it doesn't have its own 
separate repository if that is what you meant).

>  It’s code root is 
> .../swig/python/gdal-utils/osgeo_utils<https://github.com/OSGeo/gdal/tree/master/swig/python/gdal-utils/osgeo_utils>
>  . It is nested inside the gdal-utils project. In your IDE set gdal-utils as 
> the root folder. (not sure about this. Is it better to have osgeo_utils as 
> root?)
In order to make python use the package without installing it you need to add 
its parent folder (that is  
.../swig/python/gdal-utils/<https://github.com/OSGeo/gdal/tree/master/swig/python/gdal-utils/osgeo_utils>
 ) as a source dir), which practically inserts that folder at the top of your 
`sys.path`, which is where python looks for packages when you import. So 
because your `site-packages` folder is in your `sys.path` by default, then any 
package under it is recognized.
When you insert any other folder before it then it gets a higher priority (and 
that is what happens when you mark as source dir).

> Contribute changes with Pull Requests<https://github.com/OSGeo/gdal/pulls> 
> from your fork to the main GDAL project (and use gdal-utils label?)
Right, A dedicated label sounds like a good idea.

Please feel free to make a PR to improve the above docs by editing the RST 
pages (https://github.com/OSGeo/gdal/tree/master/doc/source).
Feel free to ask further questions.

Idan.


On Tue, 1 Feb 2022 at 20:51, 
mailto:matt.wil...@yukon.ca>> wrote:
Hi,

I had to read through RFC 78: gdal-utils 
package<https://gdal.org/development/rfc/rfc78_gdal_utils_package.html> a few 
times to wrap my head around what the major components are. Here’s some writing 
of what I think I understand. Is it correct?

About GDAL Python Utilities

The GDAL python utilities are included with GDAL. If you’ve installed GDAL you 
already have them. However you may want to use a newer or older version of the 
utilities without changing GDAL. This is where gdal-utils comes in.

gdal-utils: is the Gdal Python Utilities distribution. This is what you 
install. It’s home page is https://pypi.org/project/gdal-utils/ . Install it 
with `pip install gdal-utils` and it will up or downgrade the utilities in 
place.

osgeo_utils: is the python package. This is what you use in your code after 
installing, e.g. `from osgeo_utils import ...`. If you’re not writing code, 
ignore it.

(something here, perhaps names links to the utilities that are available on 
command line)

Developers

Clone or download the gdal project: https://github.com/OSGeo/gdal/

Osgeo_utils doesn’t have a top level document yet. It’s code root is 
.../swig/python/gdal-utils/osgeo_utils<https://github.com/OSGeo/gdal/tree/master/swig/python/gdal-utils/osgeo_utils>
 . It is nested inside the gdal-utils project. In your IDE set gdal-utils as 
the root folder.

(not sure about this. Is better to have osgeo_utils as root?)

Discuss via gdal-dev mailing 
list<https://lists.osgeo.org/mailman/listinfo/gdal-dev>.

Contribute changes with Pull Requests<https://github.com/OSGeo/gdal/pulls> from 
your fork to main GDAL project (and use gdal-utils label?)


Perhaps this could serve as part of a ‘gdal-utils’ project intro page 
somewhere. I think I still have write permissions to the Trac Wiki, but is that 
doesn’t feel like the right place to me.

(BW, the Nabble archive links are broken on the RFC page.)

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca<http://yukon.ca/>
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] gdal-utils: About GDAL Python Utilities

[Matt.Wilkie]

Hi,

I had to read through RFC 78: gdal-utils 
package a few 
times to wrap my head around what the major components are. Here's some writing 
of what I think I understand. Is it correct?

About GDAL Python Utilities

The GDAL python utilities are included with GDAL. If you've installed GDAL you 
already have them. However you may want to use a newer or older version of the 
utilities without changing GDAL. This is where gdal-utils comes in.

gdal-utils: is the Gdal Python Utilities distribution. This is what you 
install. It's home page is https://pypi.org/project/gdal-utils/ . Install it 
with `pip install gdal-utils` and it will up or downgrade the utilities in 
place.

osgeo_utils: is the python package. This is what you use in your code after 
installing, e.g. `from osgeo_utils import ...`. If you're not writing code, 
ignore it.

(something here, perhaps names links to the utilities that are available on 
command line)

Developers

Clone or download the gdal project: https://github.com/OSGeo/gdal/

Osgeo_utils doesn't have a top level document yet. It's code root is 
.../swig/python/gdal-utils/osgeo_utils
 . It is nested inside the gdal-utils project. In your IDE set gdal-utils as 
the root folder.

(not sure about this. Is better to have osgeo_utils as root?)

Discuss via gdal-dev mailing 
list.

Contribute changes with Pull Requests from 
your fork to main GDAL project (and use gdal-utils label?)


Perhaps this could serve as part of a 'gdal-utils' project intro page 
somewhere. I think I still have write permissions to the Trac Wiki, but is that 
doesn't feel like the right place to me.

(BW, the Nabble archive links are broken on the RFC page.)

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Nodata is None, but still has blanks?

[Matt.Wilkie]

Tiffsplit works smoothly and easily. I appreciate adding it to my long term 
tool kit. For my current project it’s great as an exploration and learning but 
as it turns out I found a more efficient route for my current project.

The out.tif that split creates for the mask band is a 3 band RGB even though it 
looks like only the first band has any info in it. That’s not really a problem 
except that cleaning it up later adds a step (in addition to re-applying the 
georeferencing). However I discovered that I actually can use gdal_translate 
out of the box by using environment settings:

set GDAL_TIFF_INTERNAL_MASK=NO
gdal_translate sample-internal-mask.tif out.tif

…will move the internal mask to an external file,  out.tif.msk, that can then 
be manipulated (a la Gimp recipe posted earlier) or discarded as needed.

It seems that using an environment setting is required. Attempts to use -co and 
--config  “GDAL_TIFF_INTERNAL_MASK=NO” resulted in a warning or error 
respectively:

Warning 6: driver GTiff does not support creation option GDAL_TIFF_INTERNAL_MASK

FAILURE: No target dataset specified.

cheers,

-Matt

From: Even Rouault 
Sent: January 18, 2022 4:21 AM
To: Frank Warmerdam ; Matt.Wilkie 
Cc: gdal dev 
Subject: Re: [gdal-dev] Nodata is None, but still has blanks?

>  I find that if I do:
tiffcp sample-no-mask.tif,0 x0.tif

I get an "x0.tif" with just the jpeg image, and not the mask.  That may be 
helpful for you (without uncompressing the jpeg image).

tiffcp doesn't use the raw interface of libtiff, hence JPEG 
decompression will occur.

tiffsplit avoids that as it does use the raw interface:

tiffsplit sample-no-mask.tif out

will generate a outaaa.tif file with the original JPEG content (the geotif tags 
will be lost however, and will have to be reinjected 
https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/samples/gdalcopyproj.py<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=https%3a%2f%2fgithub.com%2fOSGeo%2fgdal%2fblob%2fmaster%2fswig%2fpython%2fgdal%2dutils%2fosgeo%5futils%2fsamples%2fgdalcopyproj.py=B38C0EF4-D5D9-7605-8E85-D961C3CFDB4A=c132af8ee7c9d1278d61a701569070a095ce962e-99d2fc4ea472c846533a93775fde3b12fb6fdbbe>
 for example)

Even

--

http://www.spatialys.com<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=B38C0EF4-D5D9-7605-8E85-D961C3CFDB4A=c132af8ee7c9d1278d61a701569070a095ce962e-da8823f42eeaed576a737271c27f693b63c86c21>

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] Potential bug or misunderstanding. Please help

[Matt.Wilkie]

Would an in-memory /vsi also work as an intermediate?  Expressed as windows CMD 
batch file I have handy:
gdalwarp -of vrt –t_srs epsg:3578 %1 /vsistdout | gdal_translate /vsistdin/ %2


I should add that I use this like:

project-to-local.bat input.tif  outfile.tif

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

Re: [gdal-dev] Potential bug or misunderstanding. Please help

[Matt.Wilkie]

Would an in-memory /vsi also work as an intermediate?

Expressed as windows CMD batch file I have handy:

gdalwarp -of vrt –t_srs epsg:3578 %1 /vsistdout | gdal_translate /vsistdin/ %2


-Matt

From: gdal-dev  On Behalf Of Even Rouault
Sent: January 19, 2022 9:00 AM
To: Billy Babis ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Potential bug or misunderstanding. Please help


Billy,

bandList is currently ignored in the separate=True case. This would be a 
reasonable enhancement request that you could file on 
https://github.com/OSGeo/gdal/issues/new

One of the workaround for now is to postedit the generate VRT file and 
substitute the 1 element with 
2

You could also use use gdal.Translate() to produce a VRT with the second band 
of each source file, and run gdal.BuildVRT() on those files, but that would 
produce temporary files. If you don't need to serialize the output_vrt_fname, 
but use the returned dataset object with other APIs in the same process, you 
could avoid actual file creation by using an empty name for those intermediate 
VRT files, get the returned dataset object, and pass an array with those 
dataset objects instead of lst_of_input_geotiff_files.

Even

Le 19/01/2022 à 16:51, Billy Babis a écrit :
Hi there.

Given 3 input geotiffs (each with 2 bands), I want to create 1 geotiff with 3 
bands (the 2nd band from each aforementioned geotiff). This should be quite 
simple using gdal.BuildVRT and gdal.translate. But, gdal.BuildVRT will not 
allow me to select the 2nd band from each input image. The following code 
produces the following warning:

Code:
gdal.BuildVRT(output_vrt_fname, lst_of_input_geotiff_files, separate=True, 
bandList=[2])

Warning:
 "fname.tif has 2 bands. Only the first one will be taken into account for the 
-separate case"

Of course, I do not want to only use the first band, I want to use the 2nd band 
to create the VRT. Please let me know what I'm doing wrong and how I can fix it.

Best,
Billy B

--
Billy Babis
PhD Candidate: Environmental Economics & Policy (ABD)
Computer Science B.S.
1-(908)-403-2006




___

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] Nodata is None, but still has blanks?

[Matt.Wilkie]

Never mind! Paying attention to the version numbers answers the Q: Conda is 
built from libtiff 4.x while GnuWin32 libtiff is back on version 3.x.

I’ll see if I can find a contact at libtiff.org about adding a line about the 
library and tools being available via Conda as well as the other sources.

-Matt

From: gdal-dev  On Behalf Of 
matt.wil...@yukon.ca
Sent: January 18, 2022 9:34 AM
To: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Nodata is None, but still has blanks?

A little off topic for this list but people here might know: I see Tiff Tools 
for Windows binaries[0] haven’t been updated since 2006. Conda has Windows 
libtiff packages which include the tiff tools executables.[1] Is the conda 
variant actually newer than 2006 or is it just current packaging?

[0]: http://www.libtiff.org/index.html
[1]: 
https://anaconda.org/conda-forge/libtiff/4.3.0/download/win-64/libtiff-4.3.0-h0c97f57_1.tar.bz2

-Matt

From: Matt.Wilkie
Sent: January 18, 2022 8:35 AM
To: gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
Subject: RE: [gdal-dev] Nodata is None, but still has blanks?

Thank you Frank and Even.

As is so often the case on this list, I now have answers, and an education. I 
have been peripherally aware of the tiff tools from libtiff, but not enough so 
that they came to mind when seeking to puzzle out things. I’m now updated. ;-)

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

Re: [gdal-dev] Nodata is None, but still has blanks?

[Matt.Wilkie]

A little off topic for this list but people here might know: I see Tiff Tools 
for Windows binaries[0] haven’t been updated since 2006. Conda has Windows 
libtiff packages which include the tiff tools executables.[1] Is the conda 
variant actually newer than 2006 or is it just current packaging?

[0]: http://www.libtiff.org/index.html
[1]: 
https://anaconda.org/conda-forge/libtiff/4.3.0/download/win-64/libtiff-4.3.0-h0c97f57_1.tar.bz2


-Matt

From: Matt.Wilkie
Sent: January 18, 2022 8:35 AM
To: gdal-dev@lists.osgeo.org
Subject: RE: [gdal-dev] Nodata is None, but still has blanks?

Thank you Frank and Even.

As is so often the case on this list, I now have answers, and an education. I 
have been peripherally aware of the tiff tools from libtiff, but not enough so 
that they came to mind when seeking to puzzle out things. I’m now updated. ;-)

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

Re: [gdal-dev] Nodata is None, but still has blanks?

[Matt.Wilkie]

Thank you Frank and Even.

As is so often the case on this list, I now have answers, and an education. I 
have been peripherally aware of the tiff tools from libtiff, but not enough so 
that they came to mind when seeking to puzzle out things. I’m now updated. ;-)

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

[gdal-dev] Nodata is None, but still has blanks?

[Matt.Wilkie]

I'm confused. The attached image has no Nodata, but when dropped into Qgis and 
ArcGIS large areas are still drawn transparently, and there is no 4th band for 
a mask or alpha. How is this happening?

(I'm trying to get completely rid of all Nodata and Mask info, so as to start 
over creating one of those from blank slate).

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Edit geotiff in Gimp and don't lose metadata

[Matt.Wilkie]

Yes, a long time wish indeed! In discussing this with a colleague today I 
traced back my first strong desrire for such a thing back to when I got into 
making elevation models and discovering the hell that results with improperly 
set nodata. In 1994!

Now my goal is figure out how to make an external .msk file without going 
through the heavy processing time of using nearblack. Does anyone have a 
utility parameter or python snippet to create a mask file from internal Nodata 
value cells? or an arbitrary pixel value?

-Matt

From: Nikos Alexandris 
Sent: January 13, 2022 4:17 PM
To: Matt.Wilkie 
Cc: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Edit geotiff in Gimp and don't lose metadata

WoW. This is really awesome!

Whish this existed many years ago... I was using GIMP and all that magic tools 
but without resampling or altering the size of the image, and carefully 
re-applying geo-tags with the use of listgeo and geotifcp (like mentioned in 
https://marc.info/?l=grass-user=121788119621237)  :-).



On 2022-01-13 23:19, matt.wil...@yukon.ca<mailto:matt.wil...@yukon.ca> wrote:


I discovered that starting with v2.10.24 Gimp knows about and keeps geotiff 
tags intact.[0] This means we can use tools like Magic Wand fuzzy select and a 
host of other tools to quickly fix a host of image issues that are difficult to 
address using command line tools and/or code.



It gets better: Gimp reads *.msk files without any fuss. This means drastically 
reduced memory requirements for quickly making simple edits only to the mask 
layer. So a person can:



-  take a large multiband band geotiff with a mask side car file, 
squash the 3+ bands into 1 band

-  open the original .msk file in Gimp

-  drag and drop the squashed 1 band file on top the Gimp session

o   selecting import tif as a Layer

-  Lock the imported layer, set to 50% transparent (or whatever). This 
is the Visual Guide layer.

-  Paint as thou wilt on the base mask layer (Black is Nodata, White is 
Data)

-  Delete the visual guide layer

-  Export the base mask layer

o   Select tiff File Type

•  use Deflate compression

•  keep metadata (Geotiff will be disabled, that's okay)

•  use "input-rgb-geotiff.tif.msk" as file name

-  DONE!



(this is so awesome. :)



-Matt

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org<mailto: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] Edit geotiff in Gimp and don't lose metadata

[Matt.Wilkie]

-  take a large multiband band geotiff with a mask side car file, 
squash the 3+ bands into 1 band

Squashing might be unnecessary. If you open the .msk file first the document is 
set to 1-bit colour mode, then anything dropped on top is converted to the same 
mode. In effect squashing on import. The only question is whether it’s more 
efficient to let Gimp do this or to process ahead of time with gdal_translate. 
Which is faster probably depends on local conditions.

(this is so awesome. I just have to keep saying that. :)

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

Re: [gdal-dev] Edit geotiff in Gimp and don't lose metadata

[Matt.Wilkie]

I discovered that starting with v2.10.24 Gimp knows about and keeps geotiff 
tags intact.[0] This means we can use tools like Magic Wand fuzzy select and a 
host of other tools to quickly fix a host of image issues that are difficult to 
address using command line tools and/or code.

It gets better: Gimp reads *.msk files without any fuss. This means drastically 
reduced memory requirements for quickly making simple edits only to the mask 
layer. So a person can:


-  take a large multiband band geotiff with a mask side car file, 
squash the 3+ bands into 1 band

-  open the original .msk file in Gimp

-  drag and drop the squashed 1 band file on top the Gimp session

o   selecting import tif as a Layer

-  Lock the imported layer, set to 50% transparent (or whatever). This 
is the Visual Guide layer.

-  Paint as thou wilt on the base mask layer (Black is Nodata, White is 
Data)

-  Delete the visual guide layer

-  Export the base mask layer

o   Select tiff File Type

§  use Deflate compression

§  keep metadata (Geotiff will be disabled, that’s okay)

§  use “input-rgb-geotiff.tif.msk” as file name

-  DONE!

(this is so awesome. :)

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

Re: [gdal-dev] Edit geotiff in Gimp and don't lose metadata

[Matt.Wilkie]

[5]: https://gis.stackexchange.com/a/115560

Oh look at that. Paint.net has supported geotiff tags since at least 2020. 
However it dumps a stack trace when I try and open the image I used in the Gimp 
experiment. The size perhaps.

-Matt

From: gdal-dev  On Behalf Of 
matt.wil...@yukon.ca
Sent: January 13, 2022 10:49 AM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Edit geotiff in Gimp and don't lose metadata

I am s happy today!

I discovered that starting with v2.10.24 Gimp knows about and keeps geotiff 
tags intact.[0] This means we can use tools like Magic Wand fuzzy select and a 
host of other tools to quickly fix a host of image issues that are difficult to 
address using command line tools and/or code.

The big one for me at the moment is dealing with the artifact fringes 
introduced on image edges when people use jpeg compression.[1] Nearblack is 
something I keep constantly close by, but it can’t deal with areas that are 
blocked by horizontal scanning.[2] Today I used Gimp to fix the nodata area for 
5.5 GB spot image in about 20 minutes, including time with experimentation. 
True it needed ~45 GB of memory[3] so this isn’t something I’ll be doing at 
home, but it worked! [4] And I didn’t need to fart around afterwards restoring 
georeferencing.[5]

[0]: https://www.gimp.org/news/2021/03/29/gimp-2-10-24-released/
[1]: https://i.imgur.com/LZNn5Zt.png
[2]: https://i.imgur.com/avkILmX.png
[3]: https://i.imgur.com/DVLK1OB.png
[4] https://i.imgur.com/KRfxM0G.png


Some post processing is needed to get things up to our local preferences, for 
instance Gimp’s alpha channel is 15-20% larger than a gdal external mask and it 
doesn’t know my new favourite compression method of ZSTD, but this is a huge 
step forward. I could jump for joy, and maybe after I send this message I will 
;-)

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] Edit geotiff in Gimp and don't lose metadata

[Matt.Wilkie]

I am s happy today!

I discovered that starting with v2.10.24 Gimp knows about and keeps geotiff 
tags intact.[0] This means we can use tools like Magic Wand fuzzy select and a 
host of other tools to quickly fix a host of image issues that are difficult to 
address using command line tools and/or code.

The big one for me at the moment is dealing with the artifact fringes 
introduced on image edges when people use jpeg compression.[1] Nearblack is 
something I keep constantly close by, but it can't deal with areas that are 
blocked by horizontal scanning.[2] Today I used Gimp to fix the nodata area for 
5.5 GB spot image in about 20 minutes, including time with experimentation. 
True it needed ~45 GB of memory[3] so this isn't something I'll be doing at 
home, but it worked! [4] And I didn't need to fart around afterwards restoring 
georeferencing.[5]

[0]: https://www.gimp.org/news/2021/03/29/gimp-2-10-24-released/
[1]: https://i.imgur.com/LZNn5Zt.png
[2]: https://i.imgur.com/avkILmX.png
[3]: https://i.imgur.com/DVLK1OB.png
[4] https://i.imgur.com/KRfxM0G.png
[5]: https://gis.stackexchange.com/a/115560


Some post processing is needed to get things up to our local preferences, for 
instance Gimp's alpha channel is 15-20% larger than a gdal external mask and it 
doesn't know my new favourite compression method of ZSTD, but this is a huge 
step forward. I could jump for joy, and maybe after I send this message I will 
;-)

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] How to fix path for gdal .py scripts in Anaconda-Windows?

[Matt.Wilkie]

It was something of the sorts. The paths were actually fine, but the "default 
windows app" is activated when calling a file with a .py extension. So calling 
ogrmerge.py on the terminal would activate the default app instead of the 
activated python =/

If you want to manage from the command line what windows does with the .py 
extension see ASSOC and FTYPE commands


-  https://ss64.com/nt/assoc.html

-  https://ss64.com/nt/ftype.html

-  https://ss64.com/nt/fassoc.txt


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

Re: [gdal-dev] How to fix path for gdal .py scripts in Anaconda-Windows?

[Matt.Wilkie]

Also see OSGEO4W make-bat-for-py.bat. 
Old but still works. ;-)

```
for %%g in (*.py) do (
  echo @python "%%~dpnxg" %%* > %%~ng.bat
  )
```

You should be able to  insert  settings if needed for environment variables (I 
didn’t test this particular formulation):

```
@echo off
for %%g in (*.py) do (
   @echo @python "%%~dpnxg" %%*
   if not exist %%~ng.bat (
  echo “set PYTHONHOME=C:\PROGRA~1\QGIS32~1.2\apps\Python39” > %%~ng.bat
  echo “set PYTHONUTF8=1” >> %%~ng.bat
  echo @python "%%~dpnxg" %%*  >> %%~ng.bat
   )
```


-Matt

From: gdal-dev  On Behalf Of Idan Miara
Sent: January 12, 2022 2:27 PM
To: Rafael Lima 
Cc: gdal dev 
Subject: Re: [gdal-dev] How to fix path for gdal .py scripts in 
Anaconda-Windows?

Hi,

I had the same problem, here's my solution :
https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/auxiliary/batch_creator.py

Idan

On Wed, 12 Jan 2022, 23:22 Rafael Lima, 
mailto:rapl...@hotmail.com>> wrote:
I thought about uninstalling it to lol
checking that PYTHONPATH or PYTHONHOME are not defined to an unexpected value 
will be sufficient
It was something of the sorts. The paths were actually fine, but the "default 
windows app" is activated when calling a file with a .py extension. So calling 
ogrmerge.py on the terminal would activate the default app instead of the 
activated python =/
Thanks!


De: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> em 
nome de Even Rouault 
mailto:even.roua...@spatialys.com>>
Enviado: quarta-feira, 12 de janeiro de 2022 16:09
Para: gdal-dev@lists.osgeo.org 
mailto:gdal-dev@lists.osgeo.org>>
Assunto: Re: [gdal-dev] How to fix path for gdal .py scripts in 
Anaconda-Windows?


I would have suggested to uninstall the offending proprietary software :-), but 
perhaps checking that PYTHONPATH or PYTHONHOME are not defined to an unexpected 
value will be sufficient
Le 12/01/2022 à 19:24, Rafael Lima a écrit :
## Expected behavior and actual behavior.

I installed gdal on a Windows environment using anaconda and tried using one of 
the .py scripts (e.g., ogrmerge.py). I expected the script to make use of the 
gdal installed in such an environment, but it seems that the script finds 
another gdal installation (screenshot below), thus throwing an error. My 
question seems to be related to the one in 
https://stackoverflow.com/questions/69328298/cannot-run-a-py-which-requires-a-gdal-module-in-anaconda-prompt-nor-in-command,
 apparently unsolved. I was wondering what is the correct approach to make sure 
gdal scripts find the correct gdal?

## Steps to reproduce the problem.
create a conda environment:
conda create -n test1 gdal
activate the environment
conda activate test1
call a gdal .py script
ogrmerge.py

In my case, that throws an error apparently because the script uses another 
gdal installation:
[cid:part1.66CDEAD4.863B431E@spatialys.com]

## Operating system

Windows 10 Pro - 20H2 version

## GDAL version and provenance

GDAL==3.0.2
GDAL==3.4.0

Thank you!
Rafael


___

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

Re: [gdal-dev] RPCs in GeoTIFF

[Matt.Wilkie]

Thanks for the nudge. I implemented that code as a simple command line script a 
couple of months ago but neglected to put it anywhere someone else could use 
it. Fixed that today ;-)

https://gist.github.com/maphew/29a73ee7a4517663cef5081c15799fc3

```
'''Copy RPC metdata from IN raster to OUT raster

Adapted from @user7821537
https://gis.stackexchange.com/questions/264644/transfer-rpc-metadata-from-one-geotiff-to-another
'''
import os
import sys
from osgeo import gdal

gdal.UseExceptions()

if len(sys.argv) < 3:
print(f"Usage: {sys.argv[0]} [in_file] [out_file]")
sys.exit(1)

infile = sys.argv[1]# source filename and path
output = sys.argv[2]# destination file

data_in = gdal.Open(infile, gdal.GA_ReadOnly)
data_out = gdal.Open(output, gdal.GA_Update)

# get the RPCs from the first file ...
rpcs = data_in.GetMetadata('RPC')

# ... write them to the second file
data_out.SetMetadata(rpcs ,'RPC')

# de-reference the datasets, which triggers gdal to save
data_in = None
data_out = None

```

-Matt

From: gdal-dev  On Behalf Of ni hao
Sent: January 4, 2022 3:00 PM
To: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] RPCs in GeoTIFF


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



Hi list,

https://gis.stackexchange.com/questions/264644/transfer-rpc-metadata-from-one-geotiff-to-another
The above link showed me how to do it.

thanks.


From: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> on 
behalf of ni hao mailto:ni_ha...@hotmail.com>>
Sent: January 3, 2022 2:03 PM
To: gdal-dev@lists.osgeo.org 
mailto:gdal-dev@lists.osgeo.org>>
Subject: [gdal-dev] RPCs in GeoTIFF

Hi list,

I use gdal-2.4.4.
I use gdalwarp to orthorectify RADARSAT-2 images with RPCs (in product.xml).
Now I have done some manipulation and save the result in GTiff (e.g., apply 
Sigma0 LUT).
How can I copy the RPC Metadata from product.xml into GTiff, so I can use 
gdalwarp RPC orthorectification?
What are the Python commands to call ?
Do I have to use VRT?

Thank you!

ds=gdal.Open(r"C:\RADARSAT-2\2008Feb10\product.xml")
ds.GetMetadata(domain='RPC')
{'HEIGHT_OFF': '1.400e+02', 'LINE_NUM_COEFF': 
'1.045003883698053e-03 -3.909147136274753e-05 9.670521993287788e-09 
-1.540774061983927e-08 -1.230034082454030e-01 2.994775580755049e-08 
3.800851793724760e-08 -7.321141763443381e-04 -3.070513868102579e-08 
-3.263789159357597e-08 -1.074590548649705e+00 -1.422870424040082e-07 
2.625729922759320e-07 3.515681824455183e-04 1.737795244565543e-08 
-1.191967212586014e-06 -6.679635986447045e-04 9.917030647189506e-09 
1.183063120029544e-06 9.020179141273994e-07', 'SAMP_SCALE': '1759', 'LONG_OFF': 
'-5.5862000e+00', 'SAMP_DEN_COEFF': '1.000e+00 
-3.873791825160030e-04 -1.654041698986790e-06 1.045768201431545e-08 
-7.357169048971613e-04 -5.663646053748347e-06 -2.270764362304154e-08 
-7.511863551214131e-07 -3.061984244257416e-07 1.247909289740085e-07 
1.555977062556944e-03 2.290649610770613e-06 1.138900631680699e-09 
-4.788174308312700e-07 1.913987966778219e-07 -5.724932691905833e-08 
5.283021826505342e-06 -2.470297385173699e-08 -5.008669473258436e-08 
1.486018900494917e-08', 'ERR_BIAS': '5.000e-01', 'SAMP_NUM_COEFF': 
'1.448914459955711e-03 3.107453606755890e-02 -1.362345055216650e-05 
-5.545316452513627e-08 1.177884096135199e+00 2.873729811978911e-04 
-1.989082475494487e-06 -1.04193883315e-03 8.430382868139787e-06 
-1.570423196668063e-06 -2.556805573684327e-01 -1.724821081109771e-05 
4.713416494967566e-07 -9.676787333902578e-04 -2.144746372053350e-06 
8.115347578419969e-07 1.767758474109092e-04 2.353735909723349e-07 
-6.444353128930045e-06 2.385344654290232e-07', 'LONG_SCALE': 
'1.859e-01', 'LINE_DEN_COEFF': '1.000e+00 
-9.804190622020718e-08 -2.386942243966555e-07 1.801605395896790e-08 
4.160880031382031e-04 -2.095235371483683e-08 -1.659952797934781e-08 
-9.496784414468805e-07 4.758902993568164e-08 5.402628711262573e-08 
4.472124981058622e-04 1.370099419045734e-08 -1.785667248922995e-08 
7.981478756298640e-07 -7.327071351551446e-09 2.253346171933004e-07 
1.637317283006711e-06 -5.178220338393535e-08 1.880261293686099e-06 
5.550513900110941e-06', 'SAMP_OFF': '1705', 'LAT_SCALE': 
'2.073e-01', 'LAT_OFF': '3.5905800e+01', 'LINE_SCALE': 
'4229', 'LINE_OFF': '4227', 'ERR_RAND': '1.000e-01', 
'HEIGHT_SCALE': '5.010e+02'}

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

Re: [gdal-dev] follow warp with translate still?

[Matt.Wilkie]

There has been an enhancement in gdalwarp to support writing directly to 
drivers, like COG, that supports only the CreateCopy mode. So you don't need 
the 2 step process and can directly gdalwarp -of COG (it will more or less do 
the 2 steps internally)
Good news! Thanks Even. I'm glad I took the time to enquire. :)

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

[gdal-dev] follow warp with translate still?

[Matt.Wilkie]

A few years ago it was best practice to follow a gdal warp process with 
translate in order to get better compression, and even overall speed although 
extra disk write is involved. Is this still the case? and when using COG output 
driver?  e.g.

gdalwarp infile.tif interim.tif
gdal_translate -of cog -co compress=zstd  interim.tif final.tif

thanks!

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] DIMAP format questions

[Matt.Wilkie]

Hi All,

I'm getting conflicting information about the contents of some of the SPOT 
imagery we've received. Can someone here verify that this archive[0]   contains 
ortho-rectified stereo imagery? If yes can gdal use it? how?  The conversions 
I've tried so far use the RPC info but not the GCP, if indeed the files named 
GROUND_*.XML are GCP.

In the process of working with these files I've come up with number other 
questions also. They pretty much can be rolled up into: what information is 
present but NOT handled and/or carried forward by gdal?

For instance aside from the GCP question, there is a MASK folder which contains 
GML files of Clouds, Region of Interest, Saturation, Snow. I've used ogr2ogr to 
create polygons of them but have not discovered how to georeferenced them so 
they can be used with the images.


[0] I've applied a 50 pixel radius Gaussian blur to the multispectral and 300 
radius to the panchromatic TIF files. So there should be no intellectual 
property concerns about keeping this archive as a test case. I used Photoshop 
to apply the blur. If that introduced something unwanted I can reapply using 
gdal-calc or similar if someone wants to provide a command line example. All 
other files were left as is.

thanks in advance for your thoughts and time,


Files attached to this message
FilenameSizeChecksum (SHA256)
SPOT6-sample-stereo-blurred.zip 476 MB  
0f730de1e711d1cd4dbfdf84ec1ec414bc98d9c00cc2439ab8f99b1838b3b6ce

Please click on the following link to download the attachments: 
https://sft.gov.yk.ca/message/y4NYAf5MVrDPa4L10XPOcS

This email or download link can be forwarded to anyone.

The attachments are available until: Saturday, 8 January.

Message ID: y4NYAf5MVrDPa4L10XPOcS
Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Transforming Raster Points from GeoTIFF with GCPs

[Matt.Wilkie]

and https://gdal.org/programs/gdaltransform.html

-Matt

From: gdal-dev  On Behalf Of Simon Eves
Sent: December 10, 2021 12:21 PM
To: Yann-Sebastien Tremblay-Johnston 
Cc: gdal dev 
Subject: Re: [gdal-dev] Transforming Raster Points from GeoTIFF with GCPs


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



Aha! That's the bit I hadn't found. Thank you! :)

On Fri, Dec 10, 2021 at 10:38 AM Yann-Sebastien Tremblay-Johnston 
mailto:yanns.tremb...@gmail.com>> wrote:
Simon,

See the Transformer API, in particular 
https://gdal.org/api/gdal_alg.html?highlight=transformer#_CPPv419GDALTransformerFunc
 (or gdal.Transformer in Python).

Sebastien

On Fri, Dec 10, 2021 at 10:24 AM Simon Eves 
mailto:simon.e...@omnisci.com>> wrote:
So we have a GeoTiFF of Sentinel-1 imagery which `gdalinfo` reports as having a 
GCP projection with 210 points, which we'd (perhaps naively?!) like to import 
directly, transforming to WGS84 on the fly.

Obviously a simple OGRSpatialReference/OGRCoordinateTransformation transform 
isn't going to do anything, so I assume we have to do the GCP transform 
manually.

I found the GCPCoordinateTransformation class and other supporting code in 
ogr2ogr_lib.cpp as used by gdaltransform, which seems to do what we need, but 
of course it's not in the core GDAL library so we can't "just use it".

Before embarking on lifting out that code to implement an on-the-fly 
transformer in our system, I just wanted to make sure I wasn't missing any 
functionality elsewhere that would make it easier.

It may also be more efficient to just kick off a gdalwarp to pre-convert the 
file so we can just import it with no transformation required, but that seems 
ugly unless it's gonna be WAY more efficient?

Thanks,

Simon

--
[http://www2.omnisci.com/l/298412/2018-09-18/3sqpg/298412/61753/OmniSci_Email_Header2.png]

Simon Eves

Senior Graphics Engineer, Rendering Group
100 Montgomery St (5th Floor), San Francisco, CA 94104, USA



Email: simon.e...@omnisci.com | Cell:

+1 (415) 902-1996





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


--
[http://www2.omnisci.com/l/298412/2018-09-18/3sqpg/298412/61753/OmniSci_Email_Header2.png]

Simon Eves

Senior Graphics Engineer, Rendering Group
100 Montgomery St (5th Floor), San Francisco, CA 94104, USA



Email: simon.e...@omnisci.com | Cell:

+1 (415) 902-1996





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

Re: [gdal-dev] Access window out of range error with standard utils

[Matt.Wilkie]

Thanks Even!

I confirm `gdaladdo infile.ext.msk` is a valid workaround on all the problem 
files I've encountered it so far (about 7).

-Matt

From: Even Rouault 
Sent: December 9, 2021 1:28 PM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Access window out of range error with standard utils


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**




Hi,

fixed per 
https://github.com/OSGeo/gdal/pull/4965<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=https%3a%2f%2fgithub.com%2fOSGeo%2fgdal%2fpull%2f4965=FBA1B73B-D2BC-7105-950E-C186D8127588=c132af8ee7c9d1278d61a701569070a095ce962e-a9aa93d5efc34217ade976530e4f23ea2fd423a9>

you can likely workaround the issue by running gdaladdo on the .vrt.msk prior 
to COGification

Even
Le 08/12/2021 à 18:35, matt.wil...@yukon.ca<mailto:matt.wil...@yukon.ca> a 
écrit :
Hello,

Why is gdal off by one for the access window? I'm using standard utilities, no 
custom code. This is only happening (so far) on one file out of dozens.

```
$ gdal_translate -co compress=zstd -co predictor=yes -co level=17 ^
 -co bigtiff=yes -of cog INFILE.rgb-nir.vrt OUTFILE.tif
Input file size is 31918, 34252
0.ERROR 5: Access window out of range in RasterIO().  Requested
(0,0) of size 250x268 on raster of 249x267.
```

This is one step of a process where I'm taking two 3-band jpeg-in-geotiff 
images, defining a nodata mask, and creating a final 4-band cloud geotiff. 
Preceding steps to this one are:

gdal_translate -b 1 RGB.tif  b-red.vrt
gdal_translate -b 2 RGB.tif  b-grn.vrt
gdal_translate -b 3 RGB.tif  b-blu.vrt
gdal_translate -b 1 432.tif  b-nir.vrt

gdalbuildvrt -separate INFILE.rgb-nir.vrt ^
  b-red.vrt b-grn.vrt b-blu.vrt b-nir.vrt

gdal_translate INFILE.rgb-nir.vrt xx-to-mask.tif
nearblack -setmask xx-to-mask.tif
move xx-mask.tif.msk INFILE.rgb-nir.vrt.msk
del xx-mask.tif

gdal_translate INFILE.rgb-nir.vrt OUTFILE.tif

(Edited for legibility, options used are left out)

Gdalinfo reports the size of each image are the same:

$ findstr "Size" *.gdal-info
RGB.b-blu.vrt.gdal-info:Size is 31918, 34252
RGB.b-blu.vrt.gdal-info:Pixel Size = (1.500,-1.500)
RGB.b-grn.vrt.gdal-info:Size is 31918, 34252
RGB.b-grn.vrt.gdal-info:Pixel Size = (1.500,-1.500)
RGB.b-red.vrt.gdal-info:Size is 31918, 34252
RGB.b-red.vrt.gdal-info:Pixel Size = (1.500,-1.500)
432.b-nir.vrt.gdal-info:Size is 31918, 34252
432.b-nir.vrt.gdal-info:Pixel Size = (1.500,-1.500)
INFILE.rgb-nir.vrt.gdal-info:Size is 31918, 34252
INFILE.rgb-nir.vrt.gdal-info:Pixel Size = (1.500,-1.500)
INFILE.rgb-nir.vrt.msk.gdal-info:Size is 31918, 34252

The source VRT and gdalinfo reports are attached.

GDAL 3.3.2, released 2021/09/01 via Qgis 3.22 on Windows 10.

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca<http://yukon.ca/>
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.




___

gdal-dev mailing list

gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=FBA1B73B-D2BC-7105-950E-C186D8127588=c132af8ee7c9d1278d61a701569070a095ce962e-64f6b49277165076d77667c02aefeb1e8bc6fdc8>

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] Access window out of range error with standard utils

[Matt.Wilkie]

Thanks Even,

I've narrowed it down to the mask file. When I move it out of the way 
gdal_translate completes without error. Perhaps my optimization step of using 
nearblack on only the infrared band instead of all 4 is at fault?

I agree it's hard to troubleshoot when the examples are not self-contained! I 
wish I were free to share the source images. I can share the .msk file though: 
https://paste.c-net.org/StrapsMaple

(First time using this service. I know nothing about it excepting it was the 
1st result when searching for "binary pastebin")


-Matt

From: gdal-dev  On Behalf Of Even Rouault
Sent: December 8, 2021 11:08 AM
To: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Access window out of range error with standard utils


Matt,

this is likely a rounding bug somewhere, but I couldn't reproduce with what you 
provided (with GDAL master and 3.3.2 on Linux). I had to fix a few hardcoded 
paths in the VRTs and create 
SPOT7_FiftyMileCkSixtyMileR_321_150cm_19Jun2015_utm7.tif and 
SPOT7_FiftyMileCkSixtyMileR_432_150cm_19Jun2015_utm7.tif with

gdal_create SPOT7_FiftyMileCkSixtyMileR_321_150cm_19Jun2015_utm7.tif -outsize 
31918 34252 -bands 4 -co tiled=yes -co sparse_ok=yes
gdaladdo SPOT7_FiftyMileCkSixtyMileR_321_150cm_19Jun2015_utm7.tif --config 
SPARSE_OK_OVERVIEW YES   # SPARSE_OK_OVERVIEW only taken into account with 
GDAL master

I also create a .vrt.msk file with

gdal_create SPOT7_FiftyMileCkSixtyMileR_321_150cm_19Jun2015_utm7.vrt.msk 
-outsize 31912 34252 -co sparse_ok=yes -mo INTERNAL_MASK_FLAGS_1=2 -of gtiff
gdaladdo SPOT7_FiftyMileCkSixtyMileR_321_150cm_19Jun2015_utm7.vrt.msk --config 
SPARSE_OK_OVERVIEW YES

I may have missed something. It would help if your reproducer was self 
contained.

Even
Le 08/12/2021 à 18:35, matt.wil...@yukon.ca a 
écrit :
Hello,

Why is gdal off by one for the access window? I'm using standard utilities, no 
custom code. This is only happening (so far) on one file out of dozens.

```
$ gdal_translate -co compress=zstd -co predictor=yes -co level=17 ^
 -co bigtiff=yes -of cog INFILE.rgb-nir.vrt OUTFILE.tif
Input file size is 31918, 34252
0.ERROR 5: Access window out of range in RasterIO().  Requested
(0,0) of size 250x268 on raster of 249x267.
```

This is one step of a process where I'm taking two 3-band jpeg-in-geotiff 
images, defining a nodata mask, and creating a final 4-band cloud geotiff. 
Preceding steps to this one are:

gdal_translate -b 1 RGB.tif  b-red.vrt
gdal_translate -b 2 RGB.tif  b-grn.vrt
gdal_translate -b 3 RGB.tif  b-blu.vrt
gdal_translate -b 1 432.tif  b-nir.vrt

gdalbuildvrt -separate INFILE.rgb-nir.vrt ^
  b-red.vrt b-grn.vrt b-blu.vrt b-nir.vrt

gdal_translate INFILE.rgb-nir.vrt xx-to-mask.tif
nearblack -setmask xx-to-mask.tif
move xx-mask.tif.msk INFILE.rgb-nir.vrt.msk
del xx-mask.tif

gdal_translate INFILE.rgb-nir.vrt OUTFILE.tif

(Edited for legibility, options used are left out)

Gdalinfo reports the size of each image are the same:

$ findstr "Size" *.gdal-info
RGB.b-blu.vrt.gdal-info:Size is 31918, 34252
RGB.b-blu.vrt.gdal-info:Pixel Size = (1.500,-1.500)
RGB.b-grn.vrt.gdal-info:Size is 31918, 34252
RGB.b-grn.vrt.gdal-info:Pixel Size = (1.500,-1.500)
RGB.b-red.vrt.gdal-info:Size is 31918, 34252
RGB.b-red.vrt.gdal-info:Pixel Size = (1.500,-1.500)
432.b-nir.vrt.gdal-info:Size is 31918, 34252
432.b-nir.vrt.gdal-info:Pixel Size = (1.500,-1.500)
INFILE.rgb-nir.vrt.gdal-info:Size is 31918, 34252
INFILE.rgb-nir.vrt.gdal-info:Pixel Size = (1.500,-1.500)
INFILE.rgb-nir.vrt.msk.gdal-info:Size is 31918, 34252

The source VRT and gdalinfo reports are attached.

GDAL 3.3.2, released 2021/09/01 via Qgis 3.22 on Windows 10.

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.




___

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

[gdal-dev] Access window out of range error with standard utils

[Matt.Wilkie]

Hello,

Why is gdal off by one for the access window? I'm using standard utilities, no 
custom code. This is only happening (so far) on one file out of dozens.

```
$ gdal_translate -co compress=zstd -co predictor=yes -co level=17 ^
 -co bigtiff=yes -of cog INFILE.rgb-nir.vrt OUTFILE.tif
Input file size is 31918, 34252
0.ERROR 5: Access window out of range in RasterIO().  Requested
(0,0) of size 250x268 on raster of 249x267.
```

This is one step of a process where I'm taking two 3-band jpeg-in-geotiff 
images, defining a nodata mask, and creating a final 4-band cloud geotiff. 
Preceding steps to this one are:

gdal_translate -b 1 RGB.tif  b-red.vrt
gdal_translate -b 2 RGB.tif  b-grn.vrt
gdal_translate -b 3 RGB.tif  b-blu.vrt
gdal_translate -b 1 432.tif  b-nir.vrt

gdalbuildvrt -separate INFILE.rgb-nir.vrt ^
  b-red.vrt b-grn.vrt b-blu.vrt b-nir.vrt

gdal_translate INFILE.rgb-nir.vrt xx-to-mask.tif
nearblack -setmask xx-to-mask.tif
move xx-mask.tif.msk INFILE.rgb-nir.vrt.msk
del xx-mask.tif

gdal_translate INFILE.rgb-nir.vrt OUTFILE.tif

(Edited for legibility, options used are left out)

Gdalinfo reports the size of each image are the same:

$ findstr "Size" *.gdal-info
RGB.b-blu.vrt.gdal-info:Size is 31918, 34252
RGB.b-blu.vrt.gdal-info:Pixel Size = (1.500,-1.500)
RGB.b-grn.vrt.gdal-info:Size is 31918, 34252
RGB.b-grn.vrt.gdal-info:Pixel Size = (1.500,-1.500)
RGB.b-red.vrt.gdal-info:Size is 31918, 34252
RGB.b-red.vrt.gdal-info:Pixel Size = (1.500,-1.500)
432.b-nir.vrt.gdal-info:Size is 31918, 34252
432.b-nir.vrt.gdal-info:Pixel Size = (1.500,-1.500)
INFILE.rgb-nir.vrt.gdal-info:Size is 31918, 34252
INFILE.rgb-nir.vrt.gdal-info:Pixel Size = (1.500,-1.500)
INFILE.rgb-nir.vrt.msk.gdal-info:Size is 31918, 34252

The source VRT and gdalinfo reports are attached.

GDAL 3.3.2, released 2021/09/01 via Qgis 3.22 on Windows 10.

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] georeferenced but not georeferenced?

[Matt.Wilkie]

Thank you Even, especially for highlighting a DEM should be used. I 'knew' that 
but hadn't taken the time to fully ingest the meaning. It's drastically changed 
how I'm approaching this project.

-Matt

From: Even Rouault 
Sent: November 10, 2021 2:16 PM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] georeferenced but not georeferenced?

Le 10/11/2021 à 22:03, matt.wil...@yukon.ca<mailto:matt.wil...@yukon.ca> a 
écrit :

Hi all,



I'm creating images from SPOT6/7 archive DIMAP format. The resulting images 
display in the correct position in the world when combined with other sources, 
however all the usual tools say the images have no georeferencing.



* If the images have no georeferencing how do ArcMap/Pro and Qgis know to put 
them in the right place?

Those images have RPC. You can see RPC information with gdalinfo in the "RPC 
Metadata" section

listeo/geotiffcp and a number of utilities are not RPC aware. But gdalwarp is.



* How can I convert whatever-this-is into proper georeferencing so that other 
tools know how to handle them too?
You can use gdalwarp to orthorectify them.  You should generally specify -to 
RPC_DEM=some_dem_with_elevation_as_ellipsoidal_heights or -to 
RPC_HEIGHT=ellipsoidal_height. See 
https://gdal.org/api/gdal_alg.html?highlight=rpc_dem




Attached is a sample image.  It was created with (on one line):



gdal_translate  -outsize 100 0

-co compress=zstd -co predictor=yes -co level=17 -co bigtiff=yes -of cog

.\path\to\IMG_SPOT6_MS_001_A\DIM_SPOT6_MS_201609082010410_SEN_1.XML

source_tiny.tif



The biggest problem this causes me at the moment is that gdal_pansharpen.py 
output does not have any kind of goereferencing, and I can't use old faithfuls 
like listgeo, geotiffcp, or gdalcopyproj.py to copy the coordinate system 
across.



thanks in advance for your thoughts,



Matt Wilkie

Geomatics Developer & Administrator

Environment | Technology, Innovation and Mapping

T 867-667-8133 | Yukon.ca

Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.



___

gdal-dev mailing list

gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=9E426F74-D075-BC05-8F12-F338C4BD772D=c132af8ee7c9d1278d61a701569070a095ce962e-893f1a6a29bf04fdcbd6545284ebfde530cc8c97>

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] georeferenced but not georeferenced?

[Matt.Wilkie]

Hi all,

I'm creating images from SPOT6/7 archive DIMAP format. The resulting images 
display in the correct position in the world when combined with other sources, 
however all the usual tools say the images have no georeferencing. 

* If the images have no georeferencing how do ArcMap/Pro and Qgis know to put 
them in the right place?
* How can I convert whatever-this-is into proper georeferencing so that other 
tools know how to handle them too?

Attached is a sample image.  It was created with (on one line):

gdal_translate  -outsize 100 0 
-co compress=zstd -co predictor=yes -co level=17 -co bigtiff=yes -of cog 
.\path\to\IMG_SPOT6_MS_001_A\DIM_SPOT6_MS_201609082010410_SEN_1.XML 
source_tiny.tif

The biggest problem this causes me at the moment is that gdal_pansharpen.py 
output does not have any kind of goereferencing, and I can't use old faithfuls 
like listgeo, geotiffcp, or gdalcopyproj.py to copy the coordinate system 
across.

thanks in advance for your thoughts,

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Gdal-dev mail archive search sites?

[Matt.Wilkie]

> Hi Matt,
>
> mail-archive seems to be working now: 
> https://www.mail-archive.com/gdal-dev@lists.osgeo.org/

Oh that's excellent, they picked up the history as far back as April 2020 as 
well!

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

Re: [gdal-dev] Gdal-dev mail archive search sites?

[Matt.Wilkie]

> What about starting this archive?
> https://marc.info/?q=about#Add

Thanks Markus,  I’ve asked MARC about adding gdal-dev.

I’ve done the same with Mail Archive[1] but the expected subscription 
confirmation doesn’t show up[2], though there is one dating back to 2008[3]. 
Maybe the old one is getting in the way?

[1]: https://www.mail-archive.com/faq.html#newlist
[2]: https://www.mail-archive.com/archive@mail-archive.com/maillist.html
[3]: https://www.mail-archive.com/archive@mail-archive.com/msg16721.html

I will follow up with gdal-dev-owner.

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

[gdal-dev] Add a list to MARC: gdal-dev

[Matt.Wilkie]

Hi, would you please consider adding Gdal-dev to the MARC archives?

https://lists.osgeo.org/mailman/listinfo/gdal-dev

I'm not an administrator, just someone who relies on the community and is 
having difficulties consulting historical data now that Nabble has discontinued 
it's archive.

Thank you,

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Python gdal.Translate(...) and unicode in Band Descriptions.

[Matt.Wilkie]

> you could do something along:
>
> tmp_vrt = gdal.Translate("", src_ds, format = 'VRT')
> call SetColorInterpretation() / SetDescription() on tmp_vrt
> gdal.Translate(dstname, tmp_vrt, format = 'COG', ...)


Thanks Even, that works perfectly. It also revealed more depth on my on-again 
off-again troubles: In my final write I was feeding the input vsi image to 
Translate instead of the gdal dataset object I had applied the changes to. So:

gdal.Translate(outdata, VSIPATH, ...)

instead of:

gdal.Translate(outdata, DATAOBJ, ...)

The fixed and working script is attached.

-Matt


-Original Message-
From: Even Rouault  
Sent: October 14, 2021 12:11 PM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Python gdal.Translate(...) and unicode in Band 
Descriptions.

*** External email: Do not click on links or attachments except from 
trusted senders. *** 
**


you could do something along:

tmp_vrt = gdal.Translate("", src_ds, format = 'VRT')

call SetColorInterpretation() / SetDescription() on tmp_vrt

gdal.Translate(dstname, tmp_vrt, format = 'COG', ...)


Le 14/10/2021 à 20:45, matt.wil...@yukon.ca a écrit :
>> For Problem 2 there’s a new wrinkle: in spite of the warning about 
>> being unable to save the .aux.xml info into the source zip sometimes 
>> it actually does. I haven’t sorted out when it can and can’t yet.
> Solved it: the first time the script is run it saves the Color Interp and 
> Descriptions in the .aux.xml within the zip, and the output image remains 
> Undefined for these fields. In subsequent runs saving the aux fails, and the 
> output image fields for Color Interp and Description are filled out. As I've 
> been using the same input archive repeatedly while testing this sequence was 
> being masked.
>
> So now I'm at the crux I think: how to tell gdal.Translate() to use these 
> metadata without applying them to `data_in.GetRasterBand(1)` and saving first?
>
> Is there a more efficient method than:
>
> ~~~
> data_in = gdal.Open(vsipath)
> gdal.Translate(dstname, vsipath, width=w, height=h, ...)
>
> data_out = gdal.Open(dstname)
>  data_out.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
>  data_out.GetRasterBand(1).SetDescription("Red (0.625 - 0.695 microns)")
>  ...
>
> data_in = None
> data_out = None
> ~~~
>
>
> Matt
> Geomatics Developer and Administrator | Environment | T 867-667-8133 | 
> Yukon.ca
> Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.
>
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

--
https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=7AF69E2A-CE54-D505-9E07-D6028B22F810=c132af8ee7c9d1278d61a701569070a095ce962e-3d324fe99905db1c25d2b023d6680d9f015cf9ec
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] Python gdal.Translate(...) and unicode in Band Descriptions.

[Matt.Wilkie]

> For Problem 2 there’s a new wrinkle: in spite of the warning about being 
> unable to save the .aux.xml info into the source zip sometimes it actually 
> does. I haven’t sorted out when it can and can’t yet.

Solved it: the first time the script is run it saves the Color Interp and 
Descriptions in the .aux.xml within the zip, and the output image remains 
Undefined for these fields. In subsequent runs saving the aux fails, and the 
output image fields for Color Interp and Description are filled out. As I've 
been using the same input archive repeatedly while testing this sequence was 
being masked.

So now I'm at the crux I think: how to tell gdal.Translate() to use these 
metadata without applying them to `data_in.GetRasterBand(1)` and saving first?

Is there a more efficient method than:

~~~
data_in = gdal.Open(vsipath)
gdal.Translate(dstname, vsipath, width=w, height=h, ...)

data_out = gdal.Open(dstname)
data_out.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
data_out.GetRasterBand(1).SetDescription("Red (0.625 - 0.695 microns)")
...

data_in = None
data_out = None
~~~


Matt
Geomatics Developer and Administrator | Environment | T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.


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

Re: [gdal-dev] Python gdal.Translate(...) and unicode in Band Descriptions.

[Matt.Wilkie]

Re: Problem 1 - Following a thread “Unicode support in GDAL”[0] I tried adding 
“ENCODING=UTF-8” but the COG driver doesn’t support it. So I’ve taken the 
practical route out of just spelling out micron instead of using the symbol.

For Problem 2 there’s a new wrinkle: in spite of the warning about being unable 
to save the .aux.xml info into the source zip sometimes it actually does. I 
haven’t sorted out when it can and can’t yet.

[0]: https://lists.osgeo.org/pipermail/gdal-dev/2021-August/054500.html

-Matt

From: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> On 
Behalf Of matt.wil...@yukon.ca
Sent: October 13, 2021 4:39 PM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Python gdal.Translate(...) and unicode in Band Descriptions.

Hi,

The attached python script that reads from a vendor’s SPOT6 distribution 
archive and creates a cloud optimized geotiff. Reading and writing works, but 
the band description does not take properly.

Problem 1: The greek letter Mu or 
Micron
 in description is converted to something else.

In: Near-Infrared (0.760 µm – 0.890 µm)
Gdalinfo: Near-Infrared (0.760 ┬╡m ΓÇô 0.890 ┬╡m)

I’m quite sure I’ve seen geotiff descriptions using the micron character 
before. It was seeing someone else do that the gave me the idea to do the same 
with our data. Maybe they weren’t using gdal though?

Problem 2: on close it’s trying to update the source. This fails because the 
source is zipped, which is good because we want to leave the source untouched. 
However I don’t want to rely on this to keep our house in order. What is the 
proper way to do this?

Example script run:

$ spot-zip-to-preview process-list-sample.csv preview 2048
/vsizip/A:\Imagery\Work\SPOT-raw\DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790.zip/DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_MS_001_A/DIM_SPOT6_MS_202008182031259_SEN_1.XML
0.10..20...50..81..100Warning 6: driver COG does not support creation option 
PHOTOMETRIC
Warning 1: Unable to save auxiliary information in 
/vsizip/A:\Imagery\Work\SPOT-raw\DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790.zip/DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_MS_001_A/DIM_SPOT6_MS_202008182031259_SEN_1.XML.aux.xml.

The relevant part of the script:

data_in = gdal.Open(vsipath)
#print(data_in.GetMetadata())

if "_MS_" in ds:
data_in.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
data_in.GetRasterBand(1).SetDescription("Red (0.625 µm – 0.695 µm)")
data_in.GetRasterBand(2).SetColorInterpretation(gdal.GCI_GreenBand)
data_in.GetRasterBand(2).SetDescription("Green (0.530 µm – 0.590 µm)")
data_in.GetRasterBand(3).SetColorInterpretation(gdal.GCI_BlueBand)
data_in.GetRasterBand(3).SetDescription("Blue (0.455 µm – 0.525 µm)")
#data_in.GetRasterBand(4).SetColorInterpretation("Gray")
data_in.GetRasterBand(3).SetDescription("Near-Infrared (0.760 µm – 0.890 µm)")
if "_P_" in ds:
#data_in.GetRasterBand(1).SetColorInterpretation("Gray")
data_in.GetRasterBand(1).SetDescription("Panchromatic")

gdal.Translate(dstname, vsipath, width=w, height=h,
format="COG",
noData=0,
creationOptions=options,
callback=progress_cb,
callback_data='.')


Thanks in advance for your thoughts.

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Gdal-dev mail archive search sites?

[Matt.Wilkie]

Thanks Jeff!

Here's a bookmarkable link that pre-arms the search box:

https://www.google.ca/search?q=site:lists.osgeo.org/pipermail/gdal-dev/ 

-Matt


-Original Message-
From: gdal-dev  On Behalf Of Jeff McKenna
Sent: October 14, 2021 9:28 AM
To: gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] Gdal-dev mail archive search sites?

Hi Matt!

I just use Google search to query the archives:
- goto https://google.ca
- in the search box enter:site:lists.osgeo.org/pipermail/gdal-dev/ 
yoursearchterm

This may not always contain the latest archives (depending on the speed of 
Google-mothership bots), but I find it very helpful.

-jeff




--
Jeff McKenna
GatewayGeo: Developers of MS4W, MapServer Consulting and Training co-founder of 
FOSS4G http://gatewaygeo.com/


On 2021-10-14 1:03 p.m., matt.wil...@yukon.ca wrote:
> Hi All,
> 
> Now that the gdal mail archive on Nabble is gone is what options are 
> there for searching through historical threads?
> 
> Thanks!
> 
> *Matt Wilkie*
> 
> Geomatics Developer & Administrator
> 
> Environment | Technology, Innovation and Mapping
> 
> T 867-667-8133 | _Yukon.ca _
> 
> /Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away./
> 
> 
___
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] Gdal-dev mail archive search sites?

[Matt.Wilkie]

Hi All,

Now that the gdal mail archive on Nabble is gone is what options are there for 
searching through historical threads?

Thanks!

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Enabling default compression for COG ?

[Matt.Wilkie]

 +1 from me. I think it's a great idea.

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

[gdal-dev] Python gdal.Translate(...) and unicode in Band Descriptions.

[Matt.Wilkie]

Hi,

The attached python script that reads from a vendor’s SPOT6 distribution 
archive and creates a cloud optimized geotiff. Reading and writing works, but 
the band description does not take properly.

Problem 1: The greek letter Mu or 
Micron
 in description is converted to something else.

In:   Near-Infrared (0.760 µm – 0.890 µm)
Gdalinfo: Near-Infrared (0.760 ┬╡m ΓÇô 0.890 ┬╡m)

I’m quite sure I’ve seen geotiff descriptions using the micron character 
before. It was seeing someone else do that the gave me the idea to do the same 
with our data. Maybe they weren’t using gdal though?

Problem 2: on close it’s trying to update the source. This fails because the 
source is zipped, which is good because we want to leave the source untouched. 
However I don’t want to rely on this to keep our house in order. What is the 
proper way to do this?

Example script run:

$ spot-zip-to-preview process-list-sample.csv preview 2048
/vsizip/A:\Imagery\Work\SPOT-raw\DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790.zip/DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_MS_001_A/DIM_SPOT6_MS_202008182031259_SEN_1.XML
0.10..20...50..81..100Warning 6: driver COG does not support creation option 
PHOTOMETRIC
Warning 1: Unable to save auxiliary information in 
/vsizip/A:\Imagery\Work\SPOT-raw\DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790.zip/DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_MS_001_A/DIM_SPOT6_MS_202008182031259_SEN_1.XML.aux.xml.

The relevant part of the script:

data_in = gdal.Open(vsipath)
#print(data_in.GetMetadata())

if "_MS_" in ds:
data_in.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
data_in.GetRasterBand(1).SetDescription("Red (0.625 µm – 0.695 µm)")
data_in.GetRasterBand(2).SetColorInterpretation(gdal.GCI_GreenBand)
data_in.GetRasterBand(2).SetDescription("Green (0.530 µm – 0.590 µm)")
data_in.GetRasterBand(3).SetColorInterpretation(gdal.GCI_BlueBand)
data_in.GetRasterBand(3).SetDescription("Blue (0.455 µm – 0.525 µm)")
#data_in.GetRasterBand(4).SetColorInterpretation("Gray")
data_in.GetRasterBand(3).SetDescription("Near-Infrared (0.760 µm – 0.890 
µm)")
if "_P_" in ds:
#data_in.GetRasterBand(1).SetColorInterpretation("Gray")
data_in.GetRasterBand(1).SetDescription("Panchromatic")

gdal.Translate(dstname, vsipath, width=w, height=h,
format="COG",
noData=0,
creationOptions=options,
callback=progress_cb,
callback_data='.')


Thanks in advance for your thoughts.

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

<>
Driver: GTiff/GeoTIFF
Files: 
preview\DS_SPOT6_202008182031259_LM1_LM1_LM1_LM1_W137N63_01790_A_16bit.tif
Size is 2048, 2009
Metadata:
  BAND_MODE=PMS
  CLOUDCOVER_MEASURE_DESC=Region of interest mask
  CLOUDCOVER_MEASURE_NAME=Area_Of_Interest (ROI)
  CLOUDCOVER_MEASURE_TYPE=AUTOMATIC
  CLOUDCOVER_QUALITY_TABLES=SPOT
  DATASET_DELIVERY_TYPE=NETWORK
  DATASET_JOB_ID=_1
  DATASET_PRODUCER_ADDRESS=Order Desk, PO Box 21038 Westview, Lethbridge 
Alberta, T1K6X4 - Canada
  DATASET_PRODUCER_CONTACT=geomatics-sa...@planet.com
  DATASET_PRODUCER_NAME=PLANET LABS GEOMATICS CORP
  DATASET_PRODUCTION_DATE=2021-03-08T19:03:04.00Z
  DATASET_PRODUCT_INFO=SPOT
  DATASET_PRODUCT_TYPE=NA
  EPHEMERIS_ACQUISITION_ORBIT_DIRECTION=DESCENDING
  EPHEMERIS_ACQUISITION_ORBIT_NUMBER=372
  EPHEMERIS_NADIR_LAT=6.3364908090553250e+01
  EPHEMERIS_NADIR_LON=-1.4142964751935833e+02
  FACILITY_PROCESSING_CENTER=DRS
  FACILITY_SOFTWARE=IMF
  GEOMETRIC_ATTITUDES_USED=CORRECTED
  GEOMETRIC_EPHEMERIS_USED=CORRECTED
  GEOMETRIC_GEOMETRIC_PROCESSING=SENSOR
  GEOMETRIC_GROUND_SETTING=false
  GEOMETRIC_VERTICAL_SETTING=false
  IMAGING_DATE=2020-08-18
  IMAGING_TIME=20:31:25.9
  INSTRUMENT=SPOT
  INSTRUMENT_INDEX=6
  MISSION=SPOT
  MISSION_INDEX=6
  PROCESSING_LEVEL=SENSOR
  SPECTRAL_PROCESSING=MS
Image Structure Metadata:
  COMPRESSION=ZSTD
  INTERLEAVE=PIXEL
  LAYOUT=COG
RPC Metadata:
  HEIGHT_OFF=733.005
  HEIGHT_SCALE=198.5
  LAT_OFF=63.47841961
  LAT_SCALE=0.08141634
  LINE_DEN_COEFF=1 -6.31710025626589e-10 9.28731965490496e-10 
1.59192829046089e-11 -5.70050726941095e-09 9.08080140226998e-12 
1.54765000335889e-10 -2.7984199093652e-12 8.31220369759856e-09 
-8.51265131982361e-13 -2.41639854356349e-11 -5.42315554850746e-11 
1.19053085831149e-09 1.5687078530531e-14 8.86653413828049e-10 
-1.73795579431899e-09 -2.32679132842458e-14 1.5126821229747e-12 
-3.02643106305099e-11 4.84777159259441e-16
  LINE_NUM_COEFF=0.000978214554070806 0.0100406858253757 -1.01035833628811 
-0.0002992175535402 0.000495467099218108 5.23626006660564e-07 
-3.70871671103164e-05 -0.00155354829290947 

[gdal-dev] Generating data-only footprints of imagery

[Matt.Wilkie]

Hi Folks,



The standard image footprint mechanism such as that produced by gdaltindex 
draws a box around the raster but that often includes large nodata areas which 
we don't care about. Generating polygons of only the data carrying area is 
something I've often wished for but not devoted the time and attention to 
solving. This summer I'm very happy to have finally arrived at repeatable 
solution to the challenge.



I've posted the below in answer to "Creating shapefile showing footprints of 
Rasters?"
 on GIS Stack Exchange. After some more refinement I'll also post some of the 
scripts which I built as part of a larger processing chain. (No promises on how 
long that might be. I have to work it around my main duties.)



Ensure images have defined nodata. If source images do not, fix with something 
like one of the below, where 0 or 255 is the supposed-to-be nodata value:

gdal_translate ... -a_nodata 0 ... outimage.vrt

gdal_edit -a_nodata 255 somefile.tif

Create small preview images for faster processing (useful when dealing with 
multi GB satellite images). Skip this step if you want to match value pixels' 
boundary precisely. I find using the target resolution -tr parameter or 
choosing an outsize percentage that approximates 2048 rows/columns works well. 
Three approaches shown, just use one.

gdal_translate -tr 185 185 vendor_image.tif myimage.tif

gdal_translate -outsize 5% 5% vendor_image.tif myimage.tif

gdal_translate -outsize 2048 0 vendor_image.tif myimage.tif

Force pixel values to range between 1 and 255, then rescale again so that data 
= 100 and all else are nodata. Use 
vrt to not conserve storage 
and be much faster also. The reason for scaling twice is to work around the 
problem of not knowing what range the values might be when we start. There's no 
meaning to 100 other than by default it draws as a midtone grey, in contrast to 
the traditional 1 which is visually indistinguishable from black.

gdal_translate -scale -ot byte myimage.tif myimage_8bit.vrt

gdal_translate -scale 1 255 100 100 myimage_8bit.vrt myimage_data_mask.vrt

Create a polygon of the data-only area image mask for all connected regions of 
pixels in the raster sharing a common pixel value, using 8-way connectedness.:

gdal_polygonize -8 myimage_data_mask.vrt .\footprints\myimage_data.shp



Enjoy! Pink is the normal extents area as created by gdaltindex and similar 
approaches. Black outline is our new data-only polygon.

Image: https://i.stack.imgur.com/Qv1im.png 

Windows CMD shell bulk processing example

cd .\previews

md .\byte

md .\data-only

for %a in (*.tif) do gdal_translate -scale -ot byte %a .\byte\%~na.vrt

for %a in (byte\*.vrt) do gdal_translate -scale 1 255 100 100 %a 
.\data-only\%~na.vrt

for %a in (data-only\*.vrt) do call gdal_polygonize -8 %a %~na.shp

Assemble the footprints you've created by the hundred into a single overview 
index. {AUTO_NAME} populates a field with the name of the source file for that 
polygon. See ogrmerge doc page for 
other options.

ogrmerge -single ^

  -src_layer_field_name Image_Name ^

  -src_layer_field_content {AUTO_NAME} ^

  -o ..\index_footprint_data.shp ^

  .\data-only\*.shp

Image: https://i.stack.imgur.com/DbaQU.png



Thank you https://gis.stackexchange.com/a/139045/108 for rescaling data-only to 
single pixel value tip. There are other ways of reaching the same end using 
gdal_calc but I find this scale approach is easier to wrap my head around.

I think I learned about rescaling without knowing the input values ahead of 
time from Robert Simmon but have lost the reference. You can't go wrong by 
consulting his General Introduction to 
GDAL
 series however.





Hopefully this will be useful for others also,

cheers.


Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] Motion: Approve Even Rouault as a contracted GDAL maintainer

[Matt.Wilkie]

As a result of our fundraising activity and development of NumFOCUS as a 
financial conduit, it is my pleasure to put forward a motion to approve Even 
Rouault as a contracted GDAL maintainer for the year 2021-2022 beginning on 
August 1st, 2021 and ending July 31st, 2022.



Excellent news! It’s so encouraging to see a strong measure of stability and 
future proofing established.

Matt
Geomatics Developer and Administrator | Environment | T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Downscaled geotiff with overviews

[Matt.Wilkie]

Watching this thread with interest. I’ve wanted to do the same, but from python.

-Matt

From: gdal-dev  On Behalf Of Javier Jimenez 
Shaw
Sent: July 30, 2021 11:19 AM
To: gdal dev 
Subject: [gdal-dev] Downscaled geotiff with overviews

Hi

I want to make a (power-of-two) downscaled version of a geotiff with overviews 
efficiently with C++ (not from the console).

Let's say I have a geotiff with overviews (external or internal), and I want a 
copy of it, but at 50% of the size. That is just removing the "base" of the 
pyramid, but keeping all the other levels. Then voilà, there you have the 
geotiff 50% smaller, with the overviews already there, almost without any 
computation. The same for 25% (and any power of 2).

However I do not figure out how to do it efficiently in C++. If I understand 
correctly, CreateCopy is producing an image with the same size. So to downscale 
I have to create a new image, read downscaled (that will use the overviews) and 
write. But later, I have to recreate the overviews on the new image... that are 
already in the original image.

Is there any way to make it efficiently? (as we have the data there)

Thanks.
.___ ._ ..._ .. . ._.  .___ .. __ . _. . __..  ...  ._ .__
Entre dos pensamientos racionales
hay infinitos pensamientos irracionales.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Making use of SPOT imagery mask files

[Matt.Wilkie]

Aha! A small bit of progress:



$ gdaltransform -rpc 
DS_SPOT6_201308032015087_LM1_FR1_FR1_FR1_W137N65_10315\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\DIM_SPOT6_MS_201308032015087_SEN_1.XML

Enter column line values separated by space, and press Return.

0.0 57083.0

-137.766376889524 63.1915449659547 0

9652.0 0.0

-136.420879379204 66.6213864498263 0



These are the corner coordinates of the gdalinfo report, and they match the 
image as ArcMap displays it.



[cid:image002.jpg@01D78490.652687E0]


So in python we could loop through the vertices of each GML polygon coordinate 
pair and transform using the above method, but it seems there should or could 
be a more straightforward path.

-Matt


From: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> On 
Behalf Of matt.wil...@yukon.ca
Sent: July 28, 2021 2:58 PM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Making use of SPOT imagery mask files


I'm working with SPOT 6 and 7 imagery delivered in DIMAP format. I've figured 
out how to extract the multispectral and pan-chromatic bands into geotiff with 
gdal-translate so they're easier to work with.

((cross posted to 
https://gis.stackexchange.com/questions/405077/making-use-of-spot-imagery-mask-files))

In each archive there are a series of mask files in 
gml format, which I can 
read with ogrinfo and Qgis. However the mask files don't have a coordinate 
system so I can't use them with the images.

From the ogrinfo report it appears the GML are using image row and column pixel 
dimensions. (The matching source image is 9652 x 57083.)

$ ogrinfo \SPOT6_sample_roi.gml maskfeature

INFO: Open of `\SPOT6_sample_roi.gml'

  using driver `GML' successful.

Metadata:

  NAME=Area of interest mask for product id 
SPOT6_MS_201308032015087_SEN_SPOT6_20160316_1601281mdxzlrvssw12_1



Layer name: MaskFeature

Geometry: Polygon

Feature Count: 1

Extent: (1.00, 1.00) - (9653.00, 57084.00)

Layer SRS WKT:

(unknown)

gml_id: String (0.0) NOT NULL

maskType: String (18.0)

OGRFeature(MaskFeature):0

  gml_id (String) = REGION_OF_INTEREST-0

  maskType (String) = REGION_OF_INTEREST

  POLYGON ((9645.1767578125 6.41328716278076,9645.162109375 
5.32240867614746,9645.03125 4.30024194717407,9644.7841796875 
3.2741334438324,9644.5390625 2.39344930648804,9644.181640625 
1.72693908214569,9643.7666015625 1.20388793945312,9643.4931640625 
1.0,8939.99609375 1.0,1.0 4.3671669960022,1.0 28542.5,1 57084,9653 57084,9653.0 
45960.265625,9645.1767578125 6.41328716278076))

The gdalinfo report for the source 
image also 
shows pixel coordinates for the coordinate system:

Corner Coordinates:

Upper Left  (0.0,0.0)

Lower Left  (0.0,57083.0)

Upper Right ( 9652.0,0.0)

Lower Right ( 9652.0,57083.0)

Center  ( 4826.0,28541.5)

However it also has RPC metadata that seems to be enough to have Qgis and 
ArcMap/Pro display it in the right geographical location:

RPC Metadata:

  HEIGHT_OFF=500.0

  HEIGHT_SCALE=500.0

  LAT_OFF=64.90742355

  LAT_SCALE=1.71583845

  ...snip...

  SAMP_OFF=4825

  SAMP_SCALE=4826.0

Sample files at 
https://drive.google.com/drive/folders/119QEECJ42FKt0A9mq55rRmGsT9-5nSZq?usp=sharing
 (The image has been resized to 10% of it's original size.)

How might I marry the raster coordinate system info to the mask files so I can 
use them together?
Thanks in advance for your time and thoughts,

Matt Wilkie
A / Manager
(Geomatics Developer & Administrator)
Environment | Technology, Innovation and Mapping
T 867-667-8133 | 
Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] Making use of SPOT imagery mask files

[Matt.Wilkie]

I'm working with SPOT 6 and 7 imagery delivered in DIMAP format. I've figured 
out how to extract the multispectral and pan-chromatic bands into geotiff with 
gdal-translate so they're easier to work with.

((cross posted to 
https://gis.stackexchange.com/questions/405077/making-use-of-spot-imagery-mask-files))

In each archive there are a series of mask files in 
gml format, which I can 
read with ogrinfo and Qgis. However the mask files don't have a coordinate 
system so I can't use them with the images.

>From the ogrinfo report it appears the GML are using image row and column 
>pixel dimensions. (The matching source image is 9652 x 57083.)

$ ogrinfo \SPOT6_sample_roi.gml maskfeature

INFO: Open of `\SPOT6_sample_roi.gml'

  using driver `GML' successful.

Metadata:

  NAME=Area of interest mask for product id 
SPOT6_MS_201308032015087_SEN_SPOT6_20160316_1601281mdxzlrvssw12_1



Layer name: MaskFeature

Geometry: Polygon

Feature Count: 1

Extent: (1.00, 1.00) - (9653.00, 57084.00)

Layer SRS WKT:

(unknown)

gml_id: String (0.0) NOT NULL

maskType: String (18.0)

OGRFeature(MaskFeature):0

  gml_id (String) = REGION_OF_INTEREST-0

  maskType (String) = REGION_OF_INTEREST

  POLYGON ((9645.1767578125 6.41328716278076,9645.162109375 
5.32240867614746,9645.03125 4.30024194717407,9644.7841796875 
3.2741334438324,9644.5390625 2.39344930648804,9644.181640625 
1.72693908214569,9643.7666015625 1.20388793945312,9643.4931640625 
1.0,8939.99609375 1.0,1.0 4.3671669960022,1.0 28542.5,1 57084,9653 57084,9653.0 
45960.265625,9645.1767578125 6.41328716278076))

The gdalinfo report for the source 
image also 
shows pixel coordinates for the coordinate system:

Corner Coordinates:

Upper Left  (0.0,0.0)

Lower Left  (0.0,57083.0)

Upper Right ( 9652.0,0.0)

Lower Right ( 9652.0,57083.0)

Center  ( 4826.0,28541.5)

However it also has RPC metadata that seems to be enough to have Qgis and 
ArcMap/Pro display it in the right geographical location:

RPC Metadata:

  HEIGHT_OFF=500.0

  HEIGHT_SCALE=500.0

  LAT_OFF=64.90742355

  LAT_SCALE=1.71583845

  ...snip...

  SAMP_OFF=4825

  SAMP_SCALE=4826.0

Sample files at 
https://drive.google.com/drive/folders/119QEECJ42FKt0A9mq55rRmGsT9-5nSZq?usp=sharing
 (The image has been resized to 10% of it's original size.)

How might I marry the raster coordinate system info to the mask files so I can 
use them together?
Thanks in advance for your time and thoughts,

Matt Wilkie
A / Manager
(Geomatics Developer & Administrator)
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] What is the gdal_translate scale exponent option for?

[Matt.Wilkie]

Hi Folks,

gdal_translate has the option -exponent 

 to use in conjunction with -scale that will

apply non-linear scaling with a power function. exp_val is the exponent of the 
power function (must be positive).

This describes what the function does at the mechanical level but not why or 
when it might be a good idea. What are appropriate scenarios to use this 
feature? How is it useful and what does it accomplish?
Cross-posted to 
https://gis.stackexchange.com/questions/403776/what-is-the-gdal-translate-scale-exponent-option-for


Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] DIMAP driver questions

[Matt.Wilkie]

Hi,

The DIMAP driver page lists names 
of the index files to look at. I don't have any of these filenames in our SPOT6 
folder tree. Is the driver description page incomplete or perhaps I've been 
given a different DIMAP version?

METADATA.DIM (DIMAP 1),
VOL_PHR.XML (DIMAP 2) or
VOL_PNEO.XML (DIMAP 2 VHR-2020)

By pointing gdalinfo at the xml files I think I've determined that the index 
file name is of the pattern "DIM_SPOT6_MS_201508152008275_SEN_1.XML". Is this 
right?

The directory tree I have is:

```
$ tree
Folder PATH listing for volume Data
Volume serial number is C100 AEC2:FADE
D:.
└───DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499
├───LIBRARY
└───PROD_SPOT6_001
├───LIBRARY
└───VOL_SPOT6_001_A
├───IMG_SPOT6_MS_001_A
│   ├───LIBRARY
│   ├───LINEAGE
│   └───MASKS
├───IMG_SPOT6_P_001_A
│   ├───LIBRARY
│   ├───LINEAGE
│   └───MASKS
└───LIBRARY
```

Gdalinfo report and filename listings attached.


The folders named "lineage" and "mask" have names that indicate useful info 
however gdalinfo report doesn't refer to them. Is there some other way to take 
advantage of them or is that the province of full blown remote sensing software?

I'm using GDAL Utilities v3.1.4 released 2020/10/20 via Qgis 3.18 on Windows.

Thanks in advance for your thoughts,

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\DELIVERY.PDF
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\INDEX.HTM
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\LIBRARY
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\license.pdf
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\SPOT_LIST.XML
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\LIBRARY\DS_LOGO.JPG
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\LIBRARY\DS_STYLE.XSL
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\INDEX.HTM
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\LIBRARY
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\SPOT_PROD.XML
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\LIBRARY\PROD_LOGO.JPG
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\LIBRARY\PROD_STYLE.XSL
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_P_001_A
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\INDEX.HTM
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\LIBRARY
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\SPOT_VOL.XML
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\DIM_SPOT6_MS_201508152008275_SEN_1.XML
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\ICON_SPOT6_MS_201508152008275_SEN_1.JPG
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C1.TFW
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C1.TIF
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C2.TFW
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C2.TIF
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C3.TFW
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C3.TIF
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C4.TFW
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R10C4.TIF
.\DS_SPOT6_201508152008275_LM1_LM1_LM1_LM1_W138N63_06499\PROD_SPOT6_001\VOL_SPOT6_001_A\IMG_SPOT6_MS_001_A\IMG_SPOT6_MS_201508152008275_SEN_1_R1C1.TFW

Re: [gdal-dev] Install GDAL on Arcmap 10.7

[Matt.Wilkie]

Hello Sayena,

Gdal doesn’t have a graphic interface, it’s all command line utilities and 
programmer libraries. However Qgis (https://qgis.org/) is very well integrated 
with Gdal and is similar to ArcMap in a lot of ways and the members of the Qgis 
Community are generally helpful 
if the docs don’t cover your situation.

Also within ArcMap you may have access to the Data Interoperability 
extension.

-Matt

From: gdal-dev  On Behalf Of Sayena
Sent: July 9, 2021 7:32 AM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Install GDAL on Arcmap 10.7


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



Hello
I want to use Census data 2000 Urbanized area - They are in the ASCII format 
and apparently need a translator to convert them to lines/ polylines.

I searched in the internet and found out about the GDAL that can do the job. 
But just don't know how to use the graphic user interface with Arcmap 10.7

I do really appreciate your help

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

Re: [gdal-dev] Lost band name during merge

[Matt.Wilkie]

Hi Lorenzo,

Attached is a short python script to copy band info from one file to another. 
It’s not a solution but might be helpful. Usage is

gdal-copy-band-info [in_file] [band #] [out_file] [band #]


cheers,

-Matt

From: gdal-dev  On Behalf Of Lorenzo Di 
Giacomo
Sent: July 9, 2021 7:30 AM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Lost band name during merge


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



Hi, i have N netcdf OR N tiff.
I want merge them all in 1 tiff.
I do that using:
Either gdal_merge.py -separate -o 
OR gdalbuildvrt AND gdal_translate ... Both approaches works but each band has 
a name like "band_1", "band_2" ecc. the original band name is lost Is there 
a way to at least keep the band name like "original_file_name_1"... Because the 
informations abouts the bands are lost and are impossibile to know which bands 
coming from which file.
Thanks!!
''' Copy band description and color interpretation type from IN raster and BAND#
to OUT raster and BAND#'''
import os
import sys
from osgeo import gdal

gdal.UseExceptions()

if len(sys.argv) < 4:
print(f"Usage: {sys.argv[0]} [in_file] [band#] [out_file] [band #]")
sys.exit(1)

infile = sys.argv[1]# source filename and path
inband = int(sys.argv[2])   # source band number
output = sys.argv[3]# destination file
outband = int(sys.argv[4])  # destination band number

data_in = gdal.Open(infile)
band_in = data_in.GetRasterBand(inband)
data_out = gdal.Open(output, gdal.GA_Update)
band_out = data_out.GetRasterBand(outband)
band_out.SetDescription(band_in.GetDescription())
band_out.SetColorInterpretation(band_in.GetColorInterpretation())

# de-reference the datasets, which triggers gdal to save
data_in = None
data_out = None

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

Re: [gdal-dev] a flag for translate stop or abort on error?

[Matt.Wilkie]

Good to know, thanks.

For the rest of the set I'm going to unpack the archives first instead using 
/vsizip. I'm making several passes over the contents and some testing has shown 
that is slowing things down a lot. If I manage to find a reproducible set I'll 
let you know.

-Matt

From: Even Rouault 

Matt,

when you get the "TIFFReadEncodedTile() failed.", it is in a failure path of 
the geotiff driver, so the processing should stop immediately. I suspect 
something later in the stack (VRT driver?) doesn't, but without a reproducer to 
investigate what's going on exactly, hard to tell

Even
Le 08/07/2021 à 18:17, matt.wil...@yukon.ca a 
écrit :
TIFFReadEncodedTile() fai

--

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] TIFFFetchDirectory:Sanity check on directory count failed

[Matt.Wilkie]

Ah, okay. As nothing seems amiss with the translated file I'll just note the 
name and error and keep going.  Thanks for the info Even.

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

[gdal-dev] TIFFFetchDirectory:Sanity check on directory count failed

[Matt.Wilkie]

Related to my previous question today, what is the nature of 
"TIFFFetchDirectory:Sanity check on directory count failed" error?

:: gdal_translate   -co compress=zstd   -co predictor=yes   -co level=17   -co 
bigtiff=yes   -of cog   -a_nodata 0 
Finlayson9488_SP7_08Aug2017_150cm.rgb-nir.vrt 
out\Finlayson9488_SP7_08Aug2017_150cm.rgb-nir.tif
Input file size is 41114, 22580
0...10...20ERROR 1: TIFFFetchDirectory:Sanity check on directory count failed, 
this is probably not a valid IFD offset
ERROR 1: TIFFReadDirectory:Failed to read directory at offset 16879620
...30...40...50...60...70...80...90...100 - done.


The translate run finished successfully and `gdalinfo -mm ...` reports no 
errors on the output file structure. Visually the image appears intact. (Is 
scanning for minmax the best way to test for a valid file? I've been using it 
that way for years but maybe it's not best practice.)

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] a flag for translate stop or abort on error?

[Matt.Wilkie]

Hi,

Is there a flag or something to instruct gdal_translate to stop or abort when 
it encounters an error?


gdal_translate   -co compress=zstd   -co predictor=yes   -co level=17   -co 
bigtiff=yes   -of cog   -a_nodata 0 
Finlayson9488_SP6_24Aug2016_150cm.rgb-nir.vrt 
out\Finlayson9488_SP6_24Aug2016_150cm.rgb-nir.tif

Input file size is 43101, 54192
0...10..ERROR 1: LZWDecode:Corrupted LZW table at scanline 0
ERROR 1: TIFFReadEncodedTile() failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 0: TIFFReadEncodedTile() 
failed.
ERROR 1: LZWDecode:Corrupted LZW table at scanline 0
ERROR 1: TIFFReadEncodedTile() failed.
...many repeats.

The next day, no changes made to files or environment:

gdal_translate   -co compress=zstd   -co predictor=yes   -co level=17   -co 
bigtiff=yes   -of cog   -a_nodata 0 
Finlayson9488_SP6_24Aug2016_150cm.rgb-nir.vrt 
out\Finlayson9488_SP6_24Aug2016_150cm.rgb-nir.tif

Input file size is 43101, 54192
0...10...20...30...40...50... (still going)


The vrt being translated points to 4 other vrt which in turn are reading from a 
network share holding zip archives of the source imagery. (I'm taking 3 bands 
from an RGB-inzip and 1 band from a 3 band Infrared-RG-inzip and creating a 
single 4 band  RGB-I cloud geotiff.)

I suspect the problem has to do with a network timeout or something but really 
the cause doesn't matter for me here. I just want gdal_translate to stop or 
abort when something like this goes wrong. Is that possible?


Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

Re: [gdal-dev] gdaladdo bigtiff err but --config BIGTIFF was used

[Matt.Wilkie]

cogger v0.0.5 now supports .ovr external overviews. It will automatically set 
the output cog to bigtiff if required, so your use-case should now be covered.

Oh that’s great, thank you!
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] python: set CacheMax to percent, MB, GB?

[Matt.Wilkie]

Ohhh. Sorry for being dense. Sometimes the hammer has to strike a couple 
times... ;-) Thanks!
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] python: set CacheMax to percent, MB, GB?

[Matt.Wilkie]

you could also call gdal.SetConfigOption('GDAL_CACHEMAX', string), but you need 
to do that *before* reading any raster in the process (or any call to 
gdal.SetCacheMax/GetCacheMax), otherwise it will be ignored.

Thanks! ...Later: How to use it? It doesn't seem to work for me:

~~~
>>> gdal.GetCacheMax()
4294967296
>>> gdal.SetCacheMax(0)
>>> gdal.GetCacheMax()
0
>>> gdal.SetConfigOption('GDAL_CACHEMAX','50%')
>>> gdal.GetCacheMax()
0
>>> gdal.SetConfigOption('GDAL_CACHEMAX','500')
>>> gdal.GetCacheMax()
0
>>> gdal.VersionInfo()
'3020200
~~~

Python 3.8.8 (default, Feb 24 2021, 15:54:32) on Windows (64 bits).



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

[gdal-dev] python: set CacheMax to percent, MB, GB?

[Matt.Wilkie]

I see in https://gdal.org/python/osgeo.gdal-module.html#SetCacheMax we set the 
value using bytes. Is there a prebuilt method to convert units like percent, 
megabytes, gigabytes to bytes or is it up to us to build one?

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

Re: [gdal-dev] gdaladdo bigtiff err but --config BIGTIFF was used

[Matt.Wilkie]

> You can create external overviews with the -ro switch

Yeah, I have those. I was trying to create internal ones instead in order to 
use cogger ;-)

(And behind all of this, a weeping and wailing and gnashing of teeth at the 
upstream person or persons who decided that creating jpegs and throwing away 
the originals was "good enough". Oh well, my history is littered with bad calls 
too.)

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

Re: [gdal-dev] gdaladdo bigtiff err but --config BIGTIFF was used

[Matt.Wilkie]

> Your initial test.tif file needs to be a bigtiff, i.e. have been
> created with -co "BIGTIFF=YES". you can then drop the bigtiff 
> config arguments of your gdaladdo call.

Ahh, thanks Thomas, that indeed is the problem. 

Unfortunately for me it means in this case I'm stuck, as the original is 
jpeg-tiff and I don't want to pay the penalty of another decompress-recompress 
cycle to change to big tiff. However at least I know what's going on now.

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

[gdal-dev] gdaladdo bigtiff err but --config BIGTIFF was used

[Matt.Wilkie]

I am getting tiff file size exceeded error but am using bigtiff creation 
options.

~~~
$ gdaladdo -r gauss ^
--config COMPRESS_OVERVIEW WEBP ^
--config PHOTOMETRIC_OVERVIEW YCBCR ^
--config INTERLEAVE_OVERVIEW PIXEL ^
--config QUALITY 87 ^
--config BIGTIFF YES ^
--config BIGTIFF_OVERVIEW YES ^
test.tif

0...10...20...ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use 
BIGTIFF=YES creation option.
ERROR 1: An error occurred while writing a dirty block from 
GDALRasterBand::RasterIO
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES 
creation option.
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES 
creation option.
ERROR 1: WriteEncodedTile/Strip() failed.
30...40...50...60...70...80...90...100 - done.
~~~

~~~
2021-06-09  08:35 AM 4,294,967,034 test.tif
2021-03-16  09:14 AM21,765,697 test.tif.msk
~~~

Gdalinfo report attached.

GDAL 3.1.4, released 2020/10/20, Win10.

Driver: GTiff/GeoTIFF
Files: test.tif
   test.tif.msk
Size is 126015, 68149
Coordinate System is:
PROJCRS["NAD_1983_Yukon_Albers",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["Albers Equal Area",
METHOD["Albers Equal Area",
ID["EPSG",9822]],
PARAMETER["Latitude of false origin",59,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",-132.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",61.7,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",68,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",50,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",50,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",3578]]
Data axis to CRS axis mapping: 1,2
Origin = (293217.000,1126371.000)
Pixel Size = (1.500,-1.500)
Metadata:
  AREA_OR_POINT=Area
  DataType=Processed
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  293217.000, 1126371.000) (136d49'19.22"W, 64d33'30.27"N)
Lower Left  (  293217.000, 1024147.500) (136d40'52.72"W, 63d38'41.16"N)
Upper Right (  482239.500, 1126371.000) (132d52'18.43"W, 64d37'15.93"N)
Lower Right (  482239.500, 1024147.500) (132d51'34.73"W, 63d42'19.54"N)
Center  (  387728.250, 1075259.250) (134d48'36.81"W, 64d 8'43.57"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=50.198, StdDev=75.457
  NoData Value=256
  Overviews: 63008x34075, 31504x17038, 15752x8519, 7876x4260, 3938x2130, 
1969x1065, 985x533, 493x267, 247x134
  Mask Flags: PER_DATASET 
  Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=50.197691262641
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=75.456911360778
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=50.863, StdDev=75.993
  NoData Value=256
  Overviews: 63008x34075, 31504x17038, 15752x8519, 7876x4260, 3938x2130, 
1969x1065, 985x533, 493x267, 247x134
  Mask Flags: PER_DATASET 
  Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=50.862667930779
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=75.993162525726
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=44.089, StdDev=71.226
  NoData Value=256
  Overviews: 63008x34075, 31504x17038, 15752x8519, 7876x4260, 3938x2130, 
1969x1065, 985x533, 493x267, 247x134
  Mask Flags: PER_DATASET 
  Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=44.08913515331
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=71.225888610331
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] cogger output has no overviews?

[Matt.Wilkie]

Thanks Thomas.

Does cogger need a  mask band to be internal as well, and not a sidecar file?

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

Re: [gdal-dev] cogger output has no overviews?

[Matt.Wilkie]

I forgot to add: I tried adding the .ovr file in a manner similar to the Cogger 
readme example but got an error:

~~~
$ cogger -output test.tif SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif 
SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif.ovr

mucog write: load: cannot load multiple tifs if they contain overviews
~~~

-Matt

From: Matt.Wilkie
Sent: June 8, 2021 11:19 AM
To: Gdal Dev 
Subject: cogger output has no overviews?

Hi,

I took cogger for a first test drive and I'm confused. I'm relatively new to 
the world of Cloud Optimized Geotiffs, so I don't know if my confusion is about 
the COG format itself or cogger:

I fed a jpeg-in-geotiff file with external overviews (tif.ovr) and mask 
(tif.msk) to cogger. The output contains neither overviews or mask. Is that the 
expected result?

It's my understanding not having overviews is the rare case, and do that only 
if you know what you're doing (I don't!). Is there something I should be doing 
to keep or regenerate the overviews?

Gdalinfo reports for in and output attached. Gdal's validate-cog python script 
says the output is a valid cog, also attached.

Best,

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca<http://yukon.ca/>
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

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

[gdal-dev] cogger output has no overviews?

[Matt.Wilkie]

Hi,

I took cogger for a first test drive and I'm confused. I'm relatively new to 
the world of Cloud Optimized Geotiffs, so I don't know if my confusion is about 
the COG format itself or cogger:

I fed a jpeg-in-geotiff file with external overviews (tif.ovr) and mask 
(tif.msk) to cogger. The output contains neither overviews or mask. Is that the 
expected result?

It's my understanding not having overviews is the rare case, and do that only 
if you know what you're doing (I don't!). Is there something I should be doing 
to keep or regenerate the overviews?

Gdalinfo reports for in and output attached. Gdal's validate-cog python script 
says the output is a valid cog, also attached.

Best,

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Mon-Wed: Office, Thu: Remote, Fri: Away.

Driver: GTiff/GeoTIFF
Files: test.tif
Size is 126015, 68149
Coordinate System is:
PROJCRS["NAD_1983_Yukon_Albers",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["Albers Equal Area",
METHOD["Albers Equal Area",
ID["EPSG",9822]],
PARAMETER["Latitude of false origin",59,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",-132.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",61.7,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",68,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",50,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",50,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",3578]]
Data axis to CRS axis mapping: 1,2
Origin = (293217.000,1126371.000)
Pixel Size = (1.500,-1.500)
Metadata:
  AREA_OR_POINT=Area
  DataType=Processed
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  293217.000, 1126371.000) (136d49'19.22"W, 64d33'30.27"N)
Lower Left  (  293217.000, 1024147.500) (136d40'52.72"W, 63d38'41.16"N)
Upper Right (  482239.500, 1126371.000) (132d52'18.43"W, 64d37'15.93"N)
Lower Right (  482239.500, 1024147.500) (132d51'34.73"W, 63d42'19.54"N)
Center  (  387728.250, 1075259.250) (134d48'36.81"W, 64d 8'43.57"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=50.198, StdDev=75.457
  NoData Value=256
  Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=50.197691262641
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=75.456911360778
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=50.863, StdDev=75.993
  NoData Value=256
  Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=50.862667930779
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=75.993162525726
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  Min=0.000 Max=255.000 
  Minimum=0.000, Maximum=255.000, Mean=44.089, StdDev=71.226
  NoData Value=256
  Metadata:
RepresentationType=ATHEMATIC
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=44.08913515331
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=71.225888610331
Driver: GTiff/GeoTIFF
Files: SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif
   SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif.ovr
   SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif.msk.ovr
   SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif.msk
   SPOT6_321_BeaverRiver_08Sep2018_NAD83_YAlbers.tif.aux.xml
Size is 126015, 68149
Coordinate System is:
PROJCRS["NAD_1983_Yukon_Albers",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["Albers Equal Area",
METHOD["Albers Equal Area",
ID["EPSG",9822]],

Re: [gdal-dev] Introducing the cogger and godal projects

[Matt.Wilkie]

…I'd expect it to be much faster than gdal_translate -of COG (if your input is 
GeoTIFF and always properly tiled and compressed of course!).

I’m confident I know what a properly compressed geotiff is, but what is a 
properly tiled one? (or if it’s easier to define: what are improperly tiled 
geotiff?)


-Matt

From: Even Rouault 
Sent: June 4, 2021 9:33 AM
To: thomas bonfort ; Matt.Wilkie 

Cc: Gdal Dev 
Subject: Re: [gdal-dev] Introducing the cogger and godal projects


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**




I haven't benchmarked cogger but I'd expect it to be much faster than 
gdal_translate -of COG (if your input is GeoTIFF and always properly tiled and 
compressed of course!). gdal_translate -of COG uses generic GDAL API to acquire 
input data, which implies decompression / recompression. And in particular, if 
using input GeoTIFF with lossy compressed methods such as JPEG or WebP, cogger 
will avoid adding any new loss.
Le 04/06/2021 à 18:28, thomas bonfort a écrit :
I haven't extensively used -of COG (the cogger code actually predates the cog 
driver) but iirc there are at least some cases where it uses an intermediate 
file, which would imply that cogger does offer some speedups. I'll let Even 
confirm...
Regards,
Thomas
Le ven. 4 juin 2021 à 18:14, 
mailto:matt.wil...@yukon.ca>> a écrit :
Is cogger specifically for the scenario where your converting a large imagery 
library that already exists, but isn't cloud optimized? i.e. Does it offer any 
advantages over the one step `gdal_translate -of cog ...` when starting fresh?

Cheers,

Matt
Geomatics Analyst | Environment | T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.



___

gdal-dev mailing list

gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=BB1EB12E-C3F3-3E05-A4FC-DC578DED8704=c132af8ee7c9d1278d61a701569070a095ce962e-b9e494bc988145611c35e65f34295c2b143834aa>

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] Introducing the cogger and godal projects

[Matt.Wilkie]

Is cogger specifically for the scenario where your converting a large imagery 
library that already exists, but isn't cloud optimized? i.e. Does it offer any 
advantages over the one step `gdal_translate -of cog ...` when starting fresh?

Cheers,

Matt
Geomatics Analyst | Environment | T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Has gdal_ls.py been replaced?

[Matt.Wilkie]

Hi,

I only see gdal_ls.py in historical versions:
https://github.com/OSGeo/gdal/blob/e346a31f3119f3bb3f77dd7ddc8037b3049c2026/gdal/swig/python/samples/gdal_ls.py

I'm picking up on a very old suggestion[0] to use it as part of a chain to help 
gdaltindex generate footprints for images inside .zip archives. Before I go too 
far down that road though I figured I should see if this old sample is gone 
because it's been superseded or absorbed into something else.

[0]: 
https://gdal-dev.osgeo.narkive.com/o01sCeuU/gdaltindex-with-wildcard-tif-in-a-vsizip-myzip-zip

Matt Wilkie
Geomatics Developer & Administrator
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

Re: [gdal-dev] gdaldem color-relief -exact_color_entry results in RGB 000

[Matt.Wilkie]

A guess; the tiff is actually stored as 32 bit float so the values have 
decimals, and thus don’t exactly match the colour table. Maybe the Qgis layer 
is set to report without decimals?

Hmm, maybe not. Something similar was reported in GIS stack a few years ago:
https://gis.stackexchange.com/questions/197660/gdaldem-exact-color-matching-fails


Matt
Geomatics Analyst | Environment | T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

From: gdal-dev  On Behalf Of paul.m...@lfv.se
Sent: May 17, 2021 12:24 AM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] gdaldem color-relief -exact_color_entry results in RGB 000


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



Hi list, I’ve tried to use exact_color_entry on a DEM-tif file.
gdaldem color-relief 55.tif colornr.cpt 55_colrel.tif -of GTiff 
-exact_color_entry -co COMPRESS=LZW -b 1
Resulting in a blank rastefile (0 0 0).
I’ve picked (info) the colors in QGIS and there are a lot of pixels in the span 
of the color table below and QGis is presenting the values as integers.

It works fine with Nearest_color_entry.
Any ideas?
Kind regards,
Paul

colornr.cpt looks like this:
40 0 255 255
41 0 255 255
43 255 255 0
44 0 255 255
45 255 0 0
46 255 255 0
47 0 255 255
48 255 0 0
49 255 255 0
50 0 255 0
51 0 0 255
52 0 0 255
53 0 0 255
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Gdal doesn't see NaN values as NoData?

[Matt.Wilkie]

>> a) Why does gdalinfo not report the presence of nan? Meaning: how do I know 
>> when I will need to invoke (b)?
>
>It doesn't report nan as the nodata value, because the file lacks the 
>GDAL_NODATA tiff tag set at "nan". So the issue is on the producer side.

Ok, noted. I will try trace this backward and report there.


>> b) What do I need to do so that applications will use nan as nodata?
>
>   gdalwarp -dstnodata nan ...
>
>or if you want to remap nan to some other value
>
>   gdalwarp -dstnodata - ...

Wonderful, thank you.

>In most cases, you'd need to explicitly add the srcnodata value with 
>-srcnodata nan , but nan is a special case. When gdalwarp finds a pixel 
>at NaN value, be it the declared nodata value or not, it assumes this is 
>nodata, hence you can omit -srcnodata nan. Wouldn't be true if 0 or some 
>other value was the undeclared nodata  value.

Hmmm. Part of me expects gdalwarp to then record it found and used undeclared 
nodata by setting the "nodata=nan" in the metadata. Regardless thanks for this 
important detail.

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

Re: [gdal-dev] Gdal doesn't see NaN values as NoData?

[Matt.Wilkie]

Hi Folks,

I have a 32bit NDVI image that I'm told was exported from Google Earth Engine, 
with possibly a trip through R along the way. In ArcGIS desktop and Qgis there 
are "nan" values reported for the image collar, however gdalinfo doesn't say 
anything about them. If I project the image in gdalwarp the output image 
contains both nan AND nodata. This is rather confusing.

a) Why does gdalinfo not report the presence of nan? Meaning: how do I know 
when I will need to invoke (b)?
b) What do I need to do so that applications will use nan as nodata? 

Fig1: Projected output. Pink == nodata, red == 0, black == nan.


Here's the file  https://sft.gov.yk.ca/link/61czGFXOyyFsh387CxbTEF
(255 MB, link valid for 3 months. The image is free and libre.)

Gdalinfo report:

Size is 14742, 6761
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
    ELLIPSOID["WGS 84",6378137,298.257223563,
    LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
    ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
    AXIS["geodetic latitude (Lat)",north,
    ORDER[1],
    ANGLEUNIT["degree",0.0174532925199433]],
    AXIS["geodetic longitude (Lon)",east,
    ORDER[2],
    ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-137.935952550438884,61.086337635411589)
Pixel Size = (0.000269494585236,-0.000269494585236)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-137.9359526,  61.0863376) (137d56' 9.43"W, 61d 5'10.82"N)
Lower Left  (-137.9359526,  59.2642847) (137d56' 9.43"W, 59d15'51.43"N)
Upper Right (-133.9630634,  61.0863376) (133d57'47.03"W, 61d 5'10.82"N)
Lower Right (-133.9630634,  59.2642847) (133d57'47.03"W, 59d15'51.43"N)
Center  (-135.9495080,  60.1753112) (135d56'58.23"W, 60d10'31.12"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = NDVI
  Min=-0.888 Max=0.992
  Minimum=-0.888, Maximum=0.992, Mean=0.103, StdDev=0.141
  Overviews: 7371x3381, 3686x1691, 1843x846, 922x423, 461x212
  Metadata:
    STATISTICS_COVARIANCES=0.01992952634076865
    STATISTICS_MAXIMUM=0.9921875
    STATISTICS_MEAN=0.10253506184277
    STATISTICS_MINIMUM=-0.88755017518997
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=0.14117197434607


Gdalwarp command:

set _opt=-co compress=zstd -co level=17 -co num_threads=all_cpus -co 
predictor=yes -co bigtiff=yes

gdalwarp -wm 512 %_opt% -of cog -t_srs epsg:3579 -r med NDVI1984.tif 
t\NDVI1984.tif


GDAL 3.1.4, released 2020/10/20

Thanks!

Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | http://yukon.ca/
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

Re: [gdal-dev] Gdal_translate batch processing

[Matt.Wilkie]

Hi Gulnihal,

Something like this should work:

pushd C:\Users\HP\Desktop\modis
for /r %F in (*.hdf) do @echo gdal_translate -sds %F %~dpnF.tif
popd


This is a dry run to show what it will attempt, it won't actually do anything. 
Remove the "@echo" part to make it active. Double the % to use in a batch file 
(so "for /r %%F ...").

For help on what is being done here see:

 * https://ss64.com/nt/for_r.html
 * https://ss64.com/nt/syntax-args.html - what the "%~dpnF" part is doing.
 * https://ss64.com/nt/syntax-percent.html


Cheers,

-Matt

-Original Message-
From: gdal-dev  On Behalf Of gnihalk
Sent: April 29, 2021 5:36 AM
To: gdal-dev@lists.osgeo.org
Subject: [gdal-dev] Gdal_translate batch processing

*** External email: Do not click on links or attachments except from 
trusted senders. *** 
**


Hi you all. I have a huge MODIS dataset and I need to translate them to tiff 
but each data is in a different sub-folder.
C:\Users\HP\Desktop\modis\500889913\500889913\88258783\MOD10A2.2016322095024.hdf
C:\Users\HP\Desktop\modis\500889913\500889913\96113659\MOD10A2.A2016321.h20v05.006.2016340133152.hdf

.
the dataset sub-folders go like this...

and the code that im trying to run is:
for /R C:\Users\HP\Desktop\modis\500889913\500889913\88258783 %f IN
(MOD10A2.2016322095024.hdf) do gdal_translate.exe -sds -of gtiff %f 
C:\Users\HP\Desktop\modis\500889913\500889913\88258783\MOD10A2.2016322095024.tiff

I tried to make a loop for all subfolders that I have but I got an error 
massage. I would appreciate it if I can take your opinions about it. I am also 
not sure that I post this question in true title/section.
Thank you,
Gulnihal




--
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

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

[gdal-dev] Comprehensive list of values for ColorInterp

[Matt.Wilkie]

What are the useful and/or valid values for the `ColorInterp` setting of a 
raster band?
I have seen: Red, Green, Blue, Gray, Undefined.
Is this a complete list or can there be values like "Mask", "Infrared", or...?

[cid:image001.png@01D7263D.FF7F1CA0]


Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

Re: [gdal-dev] More info on gdalwarp resampling methods?

[Matt.Wilkie]

Issues opened:


· #3630 gdalwarp -multi should auto-enable GDAL_NUM_THREADS=ALL_CPUS 
<https://github.com/OSGeo/gdal/issues/3630>

· #3631 gdalwarp error "Failed to acquire WarpMutex in 
WarpRegion()"<https://github.com/OSGeo/gdal/issues/3631>

· #3632 gdalwarp mode resample much more expensive than 
needed<https://github.com/OSGeo/gdal/issues/3632>

-Matt

From: Even Rouault 
Sent: March 31, 2021 10:27 AM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] More info on gdalwarp resampling methods?


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**




Matt,

I see that in the -multi mode, the I/O threads waits for a maximum of 10 
minutes for the computation thread to have finished. And you likely hit that 
limit because of the large value of -wm you set. There's no reason we max at 10 
minutes however. Please find a but about that, including my analysis.

For better speed, I'd suggest you add --config GDAL_NUM_THREADS ALL_CPUS (or a 
given number of threads). The -multi switch only does one thread for I/O and 
one thread for computation. By adding --config GDAL_NUM_THREADS xxx you'll also 
multi-thread the computation part. The gain of -multi alone is generally 
neglectable. We should probably make -multi turn on --config GDAL_NUM_THREADS 
ALL_CPUS as well.

I also see that the mode resampling spend an awful lot of time on line 
"memset(panVals, 0, nBins*sizeof(int));" where nBins=65536 on Int16 data, and 
that for each output pixel. We could probably do something much less expensive 
for most cases. Please file a separate enhancement ticket about that particular 
case.

Even


Le 31/03/2021 à 18:46, matt.wil...@yukon.ca<mailto:matt.wil...@yukon.ca> a 
écrit :
gdalwarp --config GDAL_CACHEMAX 4096 -wm 512 -multi -co TILED=YES -co 
COMPRESS=ZSTD -t_srs EPSG:3579 -r mode -co predictor=2 
tri_mergedDEM_clipped_16bit.tif tri_mergedDEM_clipped_16bit_ytalb-rs-mode.tif

--

http://www.spatialys.com<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=682B3DEE-BED8-6F05-8639-74A30E43DC88=c132af8ee7c9d1278d61a701569070a095ce962e-453104e8853fb736505949aed5860b4e74f3ecf8>

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] More info on gdalwarp resampling methods?

[Matt.Wilkie]

It looks like my knack for stumbling into weird edge cases is still alive and 
well! Whee.

I'm using gdal from Qgis v3.16 on Windows.

I was trying the (mode,max,min,med) resampling methods and get the error 
reliably with `mode`. I thought I got it on a 2nd or 3rd method also but didn't 
repeat that this morning. I've been playing with GDAL_CACHEMAX and -wm numbers 
so that might be a factor.

Working directory is a local SSD drive and there is 32 GB physical memory, 
generally 50% or more free. Free drive space is on the low end, 10 GB on C: and 
7 GB on D:, where the working folder is.

~~~
gdalwarp ^
--config GDAL_CACHEMAX 4096 ^
-wm 512 ^
-multi ^
-co TILED=YES ^
-co COMPRESS=ZSTD ^
-t_srs EPSG:3579 ^
-r mode ^
-co predictor=2 ^
tri_mergedDEM_clipped_16bit.tif ^
tri_mergedDEM_clipped_16bit_ytalb-rs-mode.tif
~~~

Console:
~~~
$ gdalwarp --config GDAL_CACHEMAX 4096 -wm 512 -multi -co TILED=YES -co 
COMPRESS=ZSTD -t_srs EPSG:3579 -r mode -co predictor=2 
tri_mergedDEM_clipped_16bit.tif tri_mergedDEM_clipped_16bit_ytalb-rs-mode.tif
Creating output file that is 89772P x 73148L.
Processing tri_mergedDEM_clipped_16bit.tif [1/1] : 0Using internal nodata 
values (e.g. -99) for image tri_mergedDEM_clipped_16bit.tif.
Copying nodata values from source tri_mergedDEM_clipped_16bit.tif to 
destination tri_mergedDEM_clipped_16bit_ytalb-rs-mode.tif.
..ERROR 1: Failed to acquire WarpMutex in WarpRegion().
.
~~~

The info report for source image is attached, which can be downloaded from 
https://sft.gov.yk.ca/link/gV2K9F4PyH1PTJLQH8fy99  until April 14. (500mb).


-Matt

From: Even Rouault 
Sent: March 30, 2021 11:39 PM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] More info on gdalwarp resampling methods?


Matt,

(leaving out a bit your general question) The below looks is an error unrelated 
to whatever resampling method you apply. This is completely unfamiliar though. 
Do you build GDAL yourself ? Did you disable threading at build time ? Or maybe 
you're enabling threading with a huge number of threads ? Does that happen all 
the time ?

Even
Background: I'm getting error messages when trying to use some of the methods 
with 16bit imagery and I don't know if I'm seeing a gdal bug, a data bug, or 
just "you can't use that method with this raster type".

ERROR 1: Failed to acquire WarpMutex in WarpRegion().

Thanks!

Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | 
Yukon.ca<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fyukon.ca=AC1EC954-BECF-5E05-B6CB-AE995C90C7E2=c132af8ee7c9d1278d61a701569070a095ce962e-90d4f63b9d26941043d059b2a373176c1c6f200f>
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.
 ___

gdal-dev mailing list

gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fwww.spatialys.com=AC1EC954-BECF-5E05-B6CB-AE995C90C7E2=c132af8ee7c9d1278d61a701569070a095ce962e-76fdcd2f1d6f12f56ce5a5d39e4c63f2993393ad>

My software is free, but my time generally not.


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

[gdal-dev] More info on gdalwarp resampling methods?

[Matt.Wilkie]

Hi folks,

I'm seeking more information on the various gdalwarp resampling methods, the 
scenarios where they are most usefully applied, and the converse: where they 
should not be used. It would be especially helpful if there were also visuals. 
For example I know what min/max/median are when working with a list of numbers 
but my understanding of what they do when applied to rasters of various types 
is quite hazy.


Background: I'm getting error messages when trying to use some of the methods 
with 16bit imagery and I don't know if I'm seeing a gdal bug, a data bug, or 
just "you can't use that method with this raster type".

ERROR 1: Failed to acquire WarpMutex in WarpRegion().

Thanks!

Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

Re: [gdal-dev] gdalbuildvrt and band descriptions

[Matt.Wilkie]

Done!

Enable gdalbuildvrt to retain metadata such as band description
https://github.com/OSGeo/gdal/issues/3627

-Matt

From: Even Rouault 
Sent: March 26, 2021 5:09 PM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: Re: [gdal-dev] gdalbuildvrt and band descriptions

Matt,

copying over band description (and other band metadata) in the -separate case 
would make sense. Please file an enhancement ticket about that

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

[gdal-dev] gdalbuildvrt and band descriptions

[Matt.Wilkie]

Hi Folks,

`gdal_translate -b 1 infile.tif outfile.vrt` will copy the band 1 description 
to the output file.

`gdalbuild.vrt` on the other hand leaves the band descriptions behind. Is there 
a way to have vrt carry the description across also?

The project behind the question is to keep only the 4 bands of interest from 6 
bands spread across two files (Red, Green, Blue and Near infrared) and merge 
them into a single output file, and keep all their metadata. VRT is a good 
intermediary except for it drops the band descriptions. (And some other 
metadata too, but I'm less concerned about that.)

Split each band into a single vrt:

gdal_translate -b 1 rgb.tif b-red.vrt
gdal_translate -b 2 rgb.tif b-grn.vrt
gdal_translate -b 3 rgb.tif b-blu.vrt
gdal_translate -b 1 infrared.tif b-nir.vrt

Stack the 4 bands we want back together:

gdalbuildvrt -separate rgb-nir.vrt b-red.vrt b-grn.vrt b-blu.vrt b-nir.vrt


Currently I'm hand editing the last stage vrt to put the description back in, 
but it's a pain and error prone. Is there a better way?

In the end I want arcgis desktop users to see something like the image on the 
left (plus the additional NIR band) instead of the right, which is what happens 
when there's no description.

[Machine generated alternative text: WorldView-3 Blue (0.450 - 0.510) um  
WorldView-3 Green (0.510 - 0.380) um  WorldView-3 Red (0.630 - 0.690) um]

[Machine generated alternative text: @ Band_l  Band_2  Band_3]


Cheers,
Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

Re: [gdal-dev] Motion: RFC 78: gdal-utils package

[Matt.Wilkie]

Thanks, that’s exactly the info I needed to get oriented.

-Matt

From: Idan Miara 
Sent: March 26, 2021 5:29 AM
To: Matt.Wilkie 
Cc: gdal dev 
Subject: Re: [gdal-dev] Motion: RFC 78: gdal-utils package

Hi Matt,

Point noted. I've updated the Summary section:

Idan
Summary

This RFC suggests to put all the GDAL python modules (formly scripts), except 
from the GDAL core SWIG bindings, into their own distribution on pypi. The GDAL 
python sub-package osgeo.utils (introduced in GDAL 3.2) would be renamed into a 
package named osgeo_utils.

The standalone python scripts from GDAL <= 3.1 were transformed to osgeo.utils 
in GDAL 3.2. For backwards compatibility these scripts still exist and function 
as tiny wrappers around the python modules. Users of these scripts would not be 
effected from this RFC as the scripts would continue to function in GDAL 3.3 in 
the same way as in GDAL <= 3.2.

To allow maximum backwards compatibility, The osgeo package (which includes the 
GDAL core SWIG bindings) and the osgeo_utils package will continue to be 
distributed in a single sdist named gdal in pypi.

In addition, a new pure python wheel distribution named gdal-utils will be 
available in pypi under the name gdal-utils.

This will allow users who wish to upgrade the utils without upgrading the 
bindings to do so with pip install --upgrade gdal-utils (see more details in 
the following sections).


On Fri, 26 Mar 2021 at 00:11, mhw-at-yg 
mailto:matt.wil...@yukon.ca>> wrote:
Having been away from gdal-dev for some years, happily existing as an end of
the line tool only user, and only recently re-engaging with the development
community I have nothing to say on the strength and validity the RFC.

I would like to comment though that it doesn't have a high altitude overview
and say what gdal-utils is and where to look at it now. For instance is it
all the files in ./gdal/swig/python/scripts/? or gdal/swig/python/osgeo/? or
even gdal/swig/python? or somewhere else I haven't looked yet. It needn't be
a lot words (there's so much there already it's hard for a newcomer like me
to make out the structural bones). One sentence and a link would probably do
it.

I am excited about the prospect of a gdal-utils sub project, whether it
remains within the core repository or is carved off into it's own. It would
be a named place where my skills have a chance of contributing little
somethings, and I could easily ignore the much larger set of things I don't
understand in GDAL proper. ;-)

cheers,



-
-Matt
--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org<mailto: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] interpreting gdalcompare results

[Matt.Wilkie]

Hi,

How do I find out more info on what the differences are from gdalcompare? What 
does "New: n" mean?

$ gdalcompare source-image.tif  xxx-001.tif
Files differ at the binary level.
Band 1 overview count difference:
  Golden: 0
  New:5
Band 2 overview count difference:
  Golden: 0
  New:5
Band 3 overview count difference:
  Golden: 0
  New:5
Differences Found: 4


I've used raster math to subtract one image from the other and the result is 
zero in all pixels.


Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

Re: [gdal-dev] SUMMARY RE: What is lost when converting 12 bit imagery to 8 bit?

[Matt.Wilkie]

The summary also posted to the GIS Stack Exchange thread 
https://gis.stackexchange.com/questions/390315/what-is-lost-when-converting-12-bit-imagery-to-8-bit/.

Also in the credits I neglected to include @nils-erik-jørgensen.

-Matt

From: Matt.Wilkie
Sent: March 23, 2021 4:27 PM
To: 'gdal dev' 
Subject: SUMMARY RE: [gdal-dev] What is lost when converting 12 bit imagery to 
8 bit?


The position I've moved to after reading through and thinking about all the 
responses here and gdal-dev<https://lists.osgeo.org/mailman/listinfo/gdal-dev> 
(here<http://osgeo-org.1560.x6.nabble.com/gdal-dev-What-is-lost-when-converting-12-bit-imagery-to-8-bit-tt5482829.html>)
 is we need to ask for both 8 & 12 bit from the vendors.

8 bit is for visual use. It can be viewed on all displays and most programs 
without modification. Since this is 90%+ of usage scenarios having them ready 
to go is definitely beneficial to us. Additionally having the processing done 
in advance by a skilled supplier will save a lot of work -- and might even be a 
better visualization than we could produce ourselves. Producing good visuals 
across scenes systematically is not trivial.

The main drawback to automated 12 to 8 bit RGB products is that they tend to be 
very suboptimal in dark or bright areas (shadows, water, snow). This is because 
the fixed white and black point anchors needed to span the scene sets can not 
account for individual scene variation.

We should also ask imagery suppliers to describe their bit-reduction process. 
This will allow us to determine how it affects our application.

In audio terms, its a bit like cutting either the bass or the treble - what is 
best for you depends on what you want to do with the music afterwards.

  *   8 bits allows numbers from 0 to 255 (2 to the power of 8 values)
  *   12 bits allows numbers from 0 to 4095 (2 to the power of 12 values)

The simplest way conversion from 12 to 8 bits is setting the last four to zero, 
and the remaining bits shifted right four places. Its like rounding decimals to 
the nearest 100 or 1000.

What's lost is fine detail. For example the ability to discriminate between 
values of 9 and 12 - if rounded to the nearest 10 they'll both be 10 in the 
data. The same happens with binary data truncated from 12 to 8 bits.

There are other methods that could be used that preserve detail in certain 
regions of the intensity - perhaps the highest and lowest parts of the image 
can be encoded with less dynamic range, taking fewer bits, by a non-linear 
scaling of the raw data. That leaves more bits to encode the data in the middle 
of the intensity region, which might be the most important part.

12 bit is for analytic use. When using imagery for analytically (e.g. 
converting pixel values to reflectance) you probably do not want the 8 bit 
product. With the 12 bit values of 0-4095 compared to 8 bit's 0-255 there is 
opportunity to do careful scene dependent conversion in a way that best brings 
out the details available in the source data. There are a lot of methods, and 
they generally require time and patience. The challenge is sometimes called 
tone 
mapping<https://en.wikipedia.org/wiki/Tone_mapping#:%7E:text=Tone%20mapping%20is%20a%20technique,a%20more%20limited%20dynamic%20range.>.

Planet has a blog post describing how to manually convert single scene imagery 
to 8bit RGB at A Hands-On Guide to Color 
Correction<https://www.planet.com/pulse/color-correction/>. It's a good article 
as it explains the theory for how certain things are done instead of merely 
giving a recipe of steps.

To get GDAL in the mix: 
gdal-translate<https://gis.stackexchange.com/questions/tagged/gdal-translate> 
can do simple tone mapping for you with the 
-scale<https://gdal.org/programs/gdal_translate.html#cmdoption-gdal_translate-scale>
 parameter.

A special note regarding SPOT6/7: 12 bit orthorectified (level 2 processing) 
does not appear to available as a standard product, so to get 8 & 12 it might 
be necessary to purchase both Standard and Primary products. See SPOT image 
user 
guide<https://earth.esa.int/eogateway/documents/20142/37627/SPOT-6-7-imagery-user-guide.pdf>,
 particularly the descriptions of different processing levels in section 2.3.



This answer summarized from the contributions of @rs_burner, @spacedman, 
@Vince, Frank Warmerdam, and Patrick Young. Thank you!


Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca<http://yukon.ca/>
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.


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

[gdal-dev] SUMMARY RE: What is lost when converting 12 bit imagery to 8 bit?

[Matt.Wilkie]

The position I've moved to after reading through and thinking about all the 
responses here and gdal-dev 
(here)
 is we need to ask for both 8 & 12 bit from the vendors.

8 bit is for visual use. It can be viewed on all displays and most programs 
without modification. Since this is 90%+ of usage scenarios having them ready 
to go is definitely beneficial to us. Additionally having the processing done 
in advance by a skilled supplier will save a lot of work -- and might even be a 
better visualization than we could produce ourselves. Producing good visuals 
across scenes systematically is not trivial.

The main drawback to automated 12 to 8 bit RGB products is that they tend to be 
very suboptimal in dark or bright areas (shadows, water, snow). This is because 
the fixed white and black point anchors needed to span the scene sets can not 
account for individual scene variation.

We should also ask imagery suppliers to describe their bit-reduction process. 
This will allow us to determine how it affects our application.

In audio terms, its a bit like cutting either the bass or the treble - what is 
best for you depends on what you want to do with the music afterwards.

  *   8 bits allows numbers from 0 to 255 (2 to the power of 8 values)
  *   12 bits allows numbers from 0 to 4095 (2 to the power of 12 values)

The simplest way conversion from 12 to 8 bits is setting the last four to zero, 
and the remaining bits shifted right four places. Its like rounding decimals to 
the nearest 100 or 1000.

What's lost is fine detail. For example the ability to discriminate between 
values of 9 and 12 - if rounded to the nearest 10 they'll both be 10 in the 
data. The same happens with binary data truncated from 12 to 8 bits.

There are other methods that could be used that preserve detail in certain 
regions of the intensity - perhaps the highest and lowest parts of the image 
can be encoded with less dynamic range, taking fewer bits, by a non-linear 
scaling of the raw data. That leaves more bits to encode the data in the middle 
of the intensity region, which might be the most important part.

12 bit is for analytic use. When using imagery for analytically (e.g. 
converting pixel values to reflectance) you probably do not want the 8 bit 
product. With the 12 bit values of 0-4095 compared to 8 bit's 0-255 there is 
opportunity to do careful scene dependent conversion in a way that best brings 
out the details available in the source data. There are a lot of methods, and 
they generally require time and patience. The challenge is sometimes called 
tone 
mapping.

Planet has a blog post describing how to manually convert single scene imagery 
to 8bit RGB at A Hands-On Guide to Color 
Correction. It's a good article 
as it explains the theory for how certain things are done instead of merely 
giving a recipe of steps.

To get GDAL in the mix: 
gdal-translate 
can do simple tone mapping for you with the 
-scale
 parameter.

A special note regarding SPOT6/7: 12 bit orthorectified (level 2 processing) 
does not appear to available as a standard product, so to get 8 & 12 it might 
be necessary to purchase both Standard and Primary products. See SPOT image 
user 
guide,
 particularly the descriptions of different processing levels in section 2.3.



This answer summarized from the contributions of @rs_burner, @spacedman, 
@Vince, Frank Warmerdam, and Patrick Young. Thank you!


Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.


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

Re: [gdal-dev] What is lost when converting 12 bit imagery to 8 bit?

[Matt.Wilkie]

Thank you all. This is good information and helps solidify my thinking: we ask 
for both.

We want to keep getting the processed to 8bit visual imagery as we don’t have 
the capacity to that in-house for the amount of data we get. However we also 
want to have 12 bit for those occasions when analysis is primary goal, and we 
do this also.

Cheers,

-Matt

From: Frank Warmerdam 
Sent: March 17, 2021 8:29 PM
To: Patrick Young 
Cc: Matt.Wilkie ; gdal dev 
Subject: Re: [gdal-dev] What is lost when converting 12 bit imagery to 8 bit?

Patrick,

FWIW, Rob's post is on the process he uses in Photoshop to prepare images for 
various venues.  For imagery published through the platform we (Planet) do not 
use per-image white-point and black-point (or we would not have day to day, and 
scene to scene consistency).  We do apply color curves Rob prepared in our 
automated process but with "fixed" black/white point which results in an 
automated 8bit RGB product that tends to be very suboptimal in dark or bright 
areas. The imagery Planet shows in our web-explorer interface is served 
from highly compressed JPEG-in-TIFF adding an additional layer of image damage. 
:-)  While that pains me, we are keeping around nearly 3 billion scenes online 
at nearly full resolution for fast visualization so some compromises have to be 
made.

Beyond nit-picking, I think my point is:
 - given 12bit "rawer" data you have the opportunity to do careful scene 
dependent conversion to 8bit in a way that best brings out the details 
available in the source data if you have the time and patience.
 - having this process done for you in advance by a skilled supplier (perhaps 
in such a way as to maintain reasonable consistency for large area coverages) 
may actually save you a lot of work if you mostly just want to fairly generic 
visualization - and it might even be a better visualization than you would do 
yourself if you aren't going to do a lot of work.

Best regards,
Frank

On Wed, Mar 17, 2021 at 11:13 PM Patrick Young 
mailto:patrick.mckendree.yo...@gmail.com>> 
wrote:
I would guess you usually see 8bit RGB images because that is what your monitor 
can display.   What is lost is a deeper question, per channel you have to 
squeeze the original [0 - 4095] pixel value range per channel down to [0 -255], 
and there are lots of ways to do it.  The problem is sometimes called tone 
mapping<https://en.wikipedia.org/wiki/Tone_mapping#:~:text=Tone%20mapping%20is%20a%20technique,a%20more%20limited%20dynamic%20range.>.
  Planet had a nice blog post describing how they manually convert their 
imagery to 8bit RGB here<https://www.planet.com/pulse/color-correction/>.  If 
you were using the imagery for analytic things (e.g. converting pixel values to 
reflectance) you'd probably not want the 8bit product.

To get GDAL in the mix, note that gdal_translate can do simple tone mapping for 
you:  
https://gdal.org/programs/gdal_translate.html#cmdoption-gdal_translate-scale

Patrick


On Wed, Mar 17, 2021 at 3:23 PM 
mailto:matt.wil...@yukon.ca>> wrote:

SPOT 6/7 satellite imagery is captured with a dynamic range of 12 bits per 
pixel per channel (ref<https://eos.com/spot-6-and-7/>). However almost all of 
the SPOT imagery I have seen in use has been 8 bits per channel, and split into 
RGB natural colour (Bands-321) and Near-infrared-RG false colour (Bands-432). 
What information is lost in this 12 to 8 bits conversion?

I'm wondering if we should be altering our request for purchase specifications 
to deliver the full bit depth.

Although I'm referencing SPOT imagery specifically the question is general and 
really applies to any satellite or sensor system.
Cross-post: 
https://gis.stackexchange.com/questions/390315/what-is-lost-when-converting-12-bit-imagery-to-8-bit

Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | 
Yukon.ca<https://imsva91-ctp.trendmicro.com:443/wis/clicktime/v1/query?url=http%3a%2f%2fyukon.ca=4069A3E7-BDC7-3505-BA97-48548791B267=c132af8ee7c9d1278d61a701569070a095ce962e-48dacb3d842116a19074c604a613050cec3e1729>
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] What is lost when converting 12 bit imagery to 8 bit?

[Matt.Wilkie]

SPOT 6/7 satellite imagery is captured with a dynamic range of 12 bits per 
pixel per channel (ref). However almost all of 
the SPOT imagery I have seen in use has been 8 bits per channel, and split into 
RGB natural colour (Bands-321) and Near-infrared-RG false colour (Bands-432). 
What information is lost in this 12 to 8 bits conversion?

I'm wondering if we should be altering our request for purchase specifications 
to deliver the full bit depth.

Although I'm referencing SPOT imagery specifically the question is general and 
really applies to any satellite or sensor system.
Cross-post: 
https://gis.stackexchange.com/questions/390315/what-is-lost-when-converting-12-bit-imagery-to-8-bit


Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.

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

[gdal-dev] lossless jpeg crop in gdal?

[Matt.Wilkie]

I recently learned that jpegtran from the JPEG library has the ability to crop 
and save a new jpeg losslessly:

These transformations are each completely lossless and reversible. The 
transformations on the image data comprise:
· eliminate non-standard application-specific data inserted by some 
image programs,
· perform certain transformations on a file, such as:
odiscarding of colour channels (conversion to greyscale),
orotating and flipping in steps of 90 degrees,
ocropping at image block borders (every 8×8 or 16×16 pixels),
orescaling.[7]
These are lossless and reversible only regarding the image data that is kept. 
Reencoding with repeated lossy quantisation of the image data (generation loss) 
does not take place.

-- https://en.wikipedia.org/wiki/Libjpeg#jpegtran


Can GDAL take advantage of this for jpeg-in-geotiff images?

Matt Wilkie
Geomatics Analyst
Environment | Technology, Innovation and Mapping
T 867-667-8133 | Yukon.ca
Hours: 08:30-16:30, Tue-Wed: Office, Thu-Fri: Remote.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdal_retile: output isn't georeferenced?

[Matt.Wilkie]

> Is the "main" or "image" file geo-referenced ?
> It may be that the .tif.msk mask files are assumed to be the same as the 
> matching .tif file.

Ahh yes that's it. The tile main image is georeferenced and it's just the mask 
on it's own which is not. Thanks for pointing that out.

-Matt

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

Re: [gdal-dev] ZSTD compression level 17 smaller than 18+?

[Matt.Wilkie]

Oh that’s so weird. Well that puts my take-away as: just use Level=17 and 
forget the rest exist!

-Matt

From: Daniel Evans 
Sent: March 15, 2021 2:36 AM
To: Matt.Wilkie ; gdal-dev@lists.osgeo.org
Subject: RE: ZSTD compression level 17 smaller than 18+?


*** External email: Do not click on links or attachments except from trusted 
senders. ***

**



I don’t think it’s particularly unusual for any compression algorithm to have 
this sort of behaviour once it’s “maxed out”.

Running the in-built zstd benchmarks from the command line without any custom 
parameters – so entirely separate of how GDAL uses it – it shows your observed 
behaviour to an even greater extent, with most levels doing worse than level 1!

Level

Size

% Smaller than Level 1

1

3154783

0.00%

2

3130807

0.76%

3

3231415

-2.43%

4

3346308

-6.07%

5

3335646

-5.73%

6

3340284

-5.88%

7

3324298

-5.37%

8

3315170

-5.08%

9

3321085

-5.27%

10

3355473

-6.36%

11

3355932

-6.38%

12

3355857

-6.37%

13

3355401

-6.36%

14

3354678

-6.34%

15

3353801

-6.31%

16

3083731

2.25%

17

3142331

0.39%

18

3147466

0.23%

19

3146084

0.28%

20

3146548

0.26%

21

3145580

0.29%

22

3142972

0.37%


Table generated from output of `zstd -b1 -e22`, “zstd command line interface 
64-bits v1.4.7”.


Dr Daniel Evans
Software Developer

From: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> On 
Behalf Of matt.wil...@yukon.ca<mailto:matt.wil...@yukon.ca>
Sent: 15 March 2021 04:04
To: gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
Subject: [gdal-dev] ZSTD compression level 17 smaller than 18+?

So this is curious: output to Cloud optimized Geotiff with ZSTD compression to 
Level 17 is smaller than levels 18-22. Or at least that's my results when 
testing compression results with Gdal 3.1 via Qgis 3.18 on Windows 10. I 
expected compression levels to get smaller with higher levels, at expense of 
longer and longer computation time. Levels 1 to 16 do that. Even more curious 
than 17 being the smallest of the lot is that each step 18 to 22 is 
progressively larger, albeit by a miniscule amount.

Here's my results:

% smaller than biggest

Bytes

Filename

Compression

Predictor

Level

24.2%

200,968,821

cog-zstd-pyes-lev17.tif

zstd

yes

17

22.8%

204,567,155

cog-zstd-pyes-lev18.tif

zstd

yes

18

22.8%

204,594,678

cog-zstd-pyes-lev19.tif

zstd

yes

19

22.8%

204,594,678

cog-zstd-pyes-lev20.tif

zstd

yes

20

22.8%

204,626,722

cog-zstd-pyes-lev21.tif

zstd

yes

21

22.8%

204,635,985

cog-zstd-pyes-lev22.tif

zstd

yes

22

22.6%

205,026,301

cog-zstd-pyes-lev16.tif

zstd

yes

16

21.5%

208,039,012

cog-zstd-pyes-lev15.tif

zstd

yes

15

21.5%

208,094,559

cog-zstd-pyes-lev13.tif

zstd

yes

13

21.5%

208,094,560

cog-zstd-pyes-lev14.tif

zstd

yes

14

21.3%

208,714,731

cog-zstd-pyes-lev12.tif

zstd

yes

12

21.2%

208,989,301

cog-zstd-pyes-lev11.tif

zstd

yes

11

21.2%

208,989,994

cog-zstd-pyes-lev10.tif

zstd

yes

10

21.1%

209,071,549

cog-zstd-pyes-lev9.tif

zstd

yes

9

21.0%

209,453,858

cog-zstd-pyes-lev8.tif

zstd

yes

8

20.9%

209,634,360

cog-zstd-pyes-lev7.tif

zstd

yes

7

20.7%

210,182,178

cog-zstd-pyes-lev6.tif

zstd

yes

6

20.4%

210,992,074

cog-zstd-pyes-lev5.tif

zstd

yes

5

18.9%

214,894,175

cog-zstd-pyes-lev4.tif

zstd

yes

4

17.4%

218,911,762

cog-zstd-pyes-lev3.tif

lzw

2

3

16.6%

220,991,236

cog-dflt-p2.tif

deflate

2

13.0%

230,553,902

cog-zstd-pyes-lev2.tif

zstd

yes

2

11.6%

234,277,147

cog-zstd-pyes-lev1.tif

zstd

yes

1

5.7%

249,979,858

cog-dflt-p1.tif

deflate

1

0.0%

265,062,854

cog-lzw-p2.tif

lzw

2



Command line pattern:


gdal_translate ^

  -co compress=zstd ^

  -co predictor=yes ^

  -co level=3 ^

  -co num_threads=all_cpus ^

  -of cog ^

  mtn-view.tif ^

  mtn-view-cog-zstd-pyes-lev3.tif


Source image:

Name


AerialPhotos/AerialPhotos


URL


https://mapservices.gov.yk.ca/imagery/rest/services/AerialPhotos/AerialPhotos/ImageServer


Source


crs='EPSG:3578' format='JPGPNG' layer='' 
url='https://mapservices.gov.yk.ca/imagery/rest/services/AerialPhotos/AerialPhotos/ImageServer'


CRS


EPSG:3578 - NAD83 / Yukon Albers - Projected


Export extent


Upper Left, Lower Left, Upper Right, Lower Right

356750.816,  702590.026

356750.816,  699097.526

361572.583,  702590.026

361572.583,  699097.526


Pixel size


0.5m



Do others see the same thing?


​Cheers,


Matt Wilkie
Environment | Information Management & Technology | 867-667-8133  
Yukon.ca<http://yukon.ca/>

T +44 (0)1756 799919
www.jbarisk.com<http://www.jbarisk.com>

[Visit our website]<http://www.jbarisk.com> 
[http://www.jbagroup.co.uk/imgstore/JBA-Email-Sig-Icons-LINKEDIN.png]  [Follow 
us on Twitter] <https://twitter.com/jbarisk>


Find out more about us here: www.jbarisk.com<http://www.jbarisk.com/> and 
foll

Re: [gdal-dev] nearblack on multi GB images very slow

[Matt.Wilkie]

>You need to count tiles in the row, not pixels.

>
>ceil(126015/128)*3*8*128*128 = 387 118 080 = half a gigabyte.

Thank you! That made a huge difference. With GDAL_CACHEMAX=1024 it finished in 
about 40 minutes.

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

[gdal-dev] ZSTD compression level 17 smaller than 18+?

[Matt.Wilkie]

So this is curious: output to Cloud optimized Geotiff with ZSTD compression to 
Level 17 is smaller than levels 18-22. Or at least that's my results when 
testing compression results with Gdal 3.1 via Qgis 3.18 on Windows 10. I 
expected compression levels to get smaller with higher levels, at expense of 
longer and longer computation time. Levels 1 to 16 do that. Even more curious 
than 17 being the smallest of the lot is that each step 18 to 22 is 
progressively larger, albeit by a miniscule amount.

Here's my results:

% smaller than biggest  Bytes   FilenameCompression Predictor   
Level
24.2%   200,968,821 cog-zstd-pyes-lev17.tif zstdyes 17
22.8%   204,567,155 cog-zstd-pyes-lev18.tif zstdyes 18
22.8%   204,594,678
cog-zstd-pyes-lev19.tif zstdyes
19
22.8%   204,594,678 cog-zstd-pyes-lev20.tif zstdyes 20
22.8%   204,626,722 cog-zstd-pyes-lev21.tif zstdyes 21
22.8%   204,635,985 cog-zstd-pyes-lev22.tif zstdyes 22
22.6%   205,026,301 cog-zstd-pyes-lev16.tif zstdyes 16
21.5%   208,039,012
cog-zstd-pyes-lev15.tif zstdyes 15
21.5%   208,094,559 cog-zstd-pyes-lev13.tif zstdyes 13
21.5%   208,094,560 cog-zstd-pyes-lev14.tif zstdyes 14
21.3%   208,714,731 cog-zstd-pyes-lev12.tif zstdyes 12
21.2%   208,989,301 cog-zstd-pyes-lev11.tif zstdyes 11
21.2%   208,989,994 cog-zstd-pyes-lev10.tif zstdyes 10
21.1%   209,071,549 cog-zstd-pyes-lev9.tif  zstdyes 9
21.0%   209,453,858 cog-zstd-pyes-lev8.tif  zstdyes 8
20.9%   209,634,360 cog-zstd-pyes-lev7.tif  zstdyes 7
20.7%   210,182,178 cog-zstd-pyes-lev6.tif  zstdyes 6
20.4%   210,992,074 cog-zstd-pyes-lev5.tif  zstdyes 5
18.9%   214,894,175 cog-zstd-pyes-lev4.tif  zstdyes 4
17.4%   218,911,762 cog-zstd-pyes-lev3.tif  lzw 2   3
16.6%   220,991,236 cog-dflt-p2.tif deflate 2
13.0%   230,553,902 cog-zstd-pyes-lev2.tif  zstdyes
2
11.6%   234,277,147 cog-zstd-pyes-lev1.tif  zstdyes 1
5.7%249,979,858 cog-dflt-p1.tif deflate 1
0.0%265,062,854 cog-lzw-p2.tif  lzw 2



Command line pattern:


gdal_translate ^

  -co compress=zstd ^

  -co predictor=yes ^

  -co level=3 ^

  -co num_threads=all_cpus ^

  -of cog ^

  mtn-view.tif ^

  mtn-view-cog-zstd-pyes-lev3.tif


Source image:


Name


AerialPhotos/AerialPhotos


URL


https://mapservices.gov.yk.ca/imagery/rest/services/AerialPhotos/AerialPhotos/ImageServer


Source


crs='EPSG:3578' format='JPGPNG' layer='' 
url='https://mapservices.gov.yk.ca/imagery/rest/services/AerialPhotos/AerialPhotos/ImageServer'


CRS


EPSG:3578 - NAD83 / Yukon Albers - Projected


Export extent


Upper Left, Lower Left, Upper Right, Lower Right

356750.816,  702590.026

356750.816,  699097.526

361572.583,  702590.026

361572.583,  699097.526


Pixel size


0.5m



Do others see the same thing?


?Cheers,

Matt Wilkie
Environment | Information Management & Technology | 867-667-8133 ?? 
Yukon.ca
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] vrt guidance: 31 band classified raster

[Matt.Wilkie]

> It was quite a bit of work to get this far and I'm not looking forward to 
> doing this 30 more times.

On this point, discovery of the day is Easy XML 
Editor which allows editing some 
attributes as if in a spreadsheet. It's alleviated much difficulty, but is 
still finicky.

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

[gdal-dev] vrt guidance: 31 band classified raster

[Matt.Wilkie]

Hello gdal-dev, it's been a long time!  I'm happy to be digging into raster 
data building again for a change, but could use some nudges in the right 
direction(s).

A few weeks ago Nasa released Landsat-derived Annual Dominant Land Cover Across 
ABoVE Core Domain, 
1984-2014.

It's composed of 175 geotiff images, with each file containing 31 bands which 
in turn correspond to a single year in the 1984-2014 period. Each band is 8bit 
unsigned integer with values from 1 to 15 and 255 as nodata. Each integer value 
corresponds to a class such as "Evergreen Forest", "Herbaceous", "Water" and so 
on.

I've been successful in manually building a VRT file using Category element for 
the classes and ColorTable entry for a palette - but only for a single band.

My question is: how to properly apply this to all bands? Do I need duplicate 
category and colortable elements to every single VRTRasterBand element or is 
there a smart way to define it once and then refer to it like a variable?  Am I 
even approaching this the right way?

It was quite a bit of work to get this far and I'm not looking forward to doing 
this 30 more times.

Sample vrt and simplified source is attached.

Matt Wilkie
Geomatics Analyst
Environment | Information Management & Technology
T 867-667-8133 | Yukon.ca



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