Re: [gdal-dev] Non able to manage properly vitual file, warping and black edges

2024-01-22 Thread andy via gdal-dev
Hi Even
first of all, thank you very much.
I was tired yesterday. I realize this, because this morning I created a new
script starting from scratch and everything seems to be working.
Yesterday I was probably looking at the wrong output folder (I am stupid as
well as tired).

It this way it works:

gdalbuildvrt -addalpha mosaic.vrt mylist

gdalwarp -dstalpha -of vrt -overwrite -s_srs EPSG:32632 -t_srs
EPSG:7791 mosaic.vrt warp.vrt

gdal_translate -b 1 -b 2 -b 3 -b 4 -co BIGTIFF=YES -co COMPRESS=JPEG -co
PHOTOMETRIC=RGB -co ALPHA=YES -co JPEG_QUALITY=75 -co TILED=YES -co
NUM_THREADS=ALL_CPUS warp.vrt output.tif

But I need another piece of advice.
Using alpha transparency, with the 4 bands, I cannot use PHOTOMETRIC=YCBCR,
and I used PHOTOMETRIC=RGB.
This way the compression is a little less effective.

Is there anything comparable, in terms of compression effectiveness, when
using alpha transparency and JPEG compression?
Should I change the type of compression?

Thank you,

Andrea

>

-- 
___

Andrea Borruso
website: https://medium.com/tantotanto
38° 7' 48" N, 13° 21' 9" E, EPSG:4326
___

"cercare e saper riconoscere chi e cosa,
 in mezzo all’inferno, non è inferno,
e farlo durare, e dargli spazio"

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


Re: [gdal-dev] Non able to manage properly vitual file, warping and black edges

2024-01-22 Thread Rahkonen Jukka via gdal-dev
Hi,

Could you add a small image to show how the result looks like? I fear my 
imagination may be inaccurate.

-Jukka Rahkonen-

Lähettäjä: gdal-dev  Puolesta andy via 
gdal-dev
Lähetetty: maanantai 22. tammikuuta 2024 18.54
Vastaanottaja: gdal dev 
Aihe: [gdal-dev] Non able to manage properly vitual file, warping and black 
edges

Hi,
I have four tif input files, these have 3 bands RGB and are JPEG compressed.
No null data set.

My goal is to mosaic these and then warp the mosaic.

I start it with:

gdalbuildvrt mosaic.vrt file1.tif  file2.tif  file3.tif  file4.tif

Then I warp the mosaic:

gdalwarp -dstalpha -of vrt -overwrite -s_srs EPSG:32632 -t_srs EPSG:7791 
mosaic.vrt warped.vrt

Then I run

gdal_translate -b 1 -b 2 -b 3 -b 4 -co BIGTIFF=YES -co COMPRESS=DEFLATE -co 
ALPHA=YES -co TILED=YES -co NUM_THREADS=ALL_CPUS  warped.vrt  output.tif

The final tiff file has black edges.
What's wrong with my procedure? How to avoid this?

Thank you

--
___

Andrea Borruso
website: https://medium.com/tantotanto
38° 7' 48" N, 13° 21' 9" E, EPSG:4326
___

"cercare e saper riconoscere chi e cosa,
 in mezzo all'inferno, non è inferno,
e farlo durare, e dargli spazio"

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


Re: [gdal-dev] Non able to manage properly vitual file, warping and black edges

2024-01-22 Thread Even Rouault via gdal-dev
ah I missed your sources were JPEG compressed. Do they have padding 
around them? If so, you might need to use nearblack to compute a mask band


Le 22/01/2024 à 18:09, andy a écrit :

Hi Even,
and thank you, I have the same result.

Do I must remove jpeg compression in gdaladdo?

Best regards



--
___

Andrea Borruso
website: https://medium.com/tantotanto
38° 7' 48" N, 13° 21' 9" E, EPSG:4326
___

"cercare e saper riconoscere chi e cosa,
 in mezzo all’inferno, non è inferno,
e farlo durare, e dargli spazio"

Italo Calvino


--
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] Non able to manage properly vitual file, warping and black edges

2024-01-22 Thread andy via gdal-dev
Hi Even,
and thank you, I have the same result.

Do I must remove jpeg compression in gdaladdo?

Best regards




-- 
___

Andrea Borruso
website: https://medium.com/tantotanto
38° 7' 48" N, 13° 21' 9" E, EPSG:4326
___

"cercare e saper riconoscere chi e cosa,
 in mezzo all’inferno, non è inferno,
e farlo durare, e dargli spazio"

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


Re: [gdal-dev] Non able to manage properly vitual file, warping and black edges

2024-01-22 Thread Even Rouault via gdal-dev


Le 22/01/2024 à 17:53, andy via gdal-dev a écrit :

Hi,
I have four tif input files, these have 3 bands RGB and are JPEG 
compressed.

No null data set.

My goal is to mosaic these and then warp the mosaic.

I start it with:

gdalbuildvrt mosaic.vrt file1.tif file2.tif file3.tif file4.tif


Try adding a -addalpha so that the mosaic.vrt is created with an alpha 
band which equals to 0 in gaps between your source tiles. The rest of 
your procedure should be fine


--
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] Non able to manage properly vitual file, warping and black edges

2024-01-22 Thread andy via gdal-dev
Hi,
I have four tif input files, these have 3 bands RGB and are JPEG compressed.
No null data set.

My goal is to mosaic these and then warp the mosaic.

I start it with:

gdalbuildvrt mosaic.vrt file1.tif  file2.tif  file3.tif  file4.tif

Then I warp the mosaic:

gdalwarp -dstalpha -of vrt -overwrite -s_srs EPSG:32632 -t_srs EPSG:7791
mosaic.vrt warped.vrt

Then I run

gdal_translate -b 1 -b 2 -b 3 -b 4 -co BIGTIFF=YES -co COMPRESS=DEFLATE -co
ALPHA=YES -co TILED=YES -co NUM_THREADS=ALL_CPUS  warped.vrt  output.tif

The final tiff file has black edges.
What's wrong with my procedure? How to avoid this?

Thank you

-- 
___

Andrea Borruso
website: https://medium.com/tantotanto
38° 7' 48" N, 13° 21' 9" E, EPSG:4326
___

"cercare e saper riconoscere chi e cosa,
 in mezzo all’inferno, non è inferno,
e farlo durare, e dargli spazio"

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


Re: [gdal-dev] warp issue: Source and target ellipsoid do not belong to the same celestial body

2024-01-22 Thread Even Rouault via gdal-dev


Le 22/01/2024 à 17:33, Wilco K via gdal-dev a écrit :

Hi,

these 2 lines work with GDAL 2.2.3, but not with GDAL 3.5.0 anymore.

gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 
+lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" -a_ullr 0.0, 
-3649.9792, 700.000, -4414.9792 
"HDF5:\"RAD_NL25_RAC_03H_202401220800.h5\"://image1/image_data" 
"test-translate.tif"


gdalwarp -t_srs EPSG:4326 -dstnodata 65535 -of GTiff 
"test-translate.tif" "test-warp.tif"



Error:
ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do 
not belong to the same celestial body
ERROR 6: Cannot find coordinate operations from 
`PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378.14,298.183263207102,LENGTHUNIT["metre",1,ID["EPSG",9001,PRIMEM["Reference 
meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122,CONVERSION["Polar 
Stereographic (variant B)",METHOD["Polar Stereographic (variant 
B)",ID["EPSG",9829]],PARAMETER["Latitude of standard 
parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude 
of 
origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False 
easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False 
northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1]]]' 
to `EPSG:4326'



But when the +a and +b in the gdal_translate call are changed, 
gdalwarp does work:
gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 
+lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0" -a_ullr 0 
-364.11191775 70.903671186 -4415003.88199764 
"HDF5:\"RAD_NL25_RAC_03H_202401220800.h5\"://image1/image_data" 
"test-translate.tif"


What is the problem?


The value of the +a and +b parameters must be in metres, not in km. So 
you add to cheat also on the geotransform. But if datum transformations 
had to be done the fact that you transform between an ellipsoid of ~ 
6000 meters to one of ~ 6000 km wouldn't work well... It probably sort 
of worked because datum transformation was skipped, but the new 
behaviour which checks that the shape of the source and target ellipsoid 
isn't too different is definitely saner and will avoid potential 
reprojection errors


--
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] warp issue: Source and target ellipsoid do not belong to the same celestial body

2024-01-22 Thread Wilco K via gdal-dev
Hi,

these 2 lines work with GDAL 2.2.3, but not with GDAL 3.5.0 anymore.

gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 +lon_0=0 
+lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" -a_ullr 0.0, -3649.9792, 
700.000, -4414.9792 
"HDF5:\"RAD_NL25_RAC_03H_202401220800.h5\"://image1/image_data" 
"test-translate.tif"

gdalwarp -t_srs EPSG:4326 -dstnodata 65535 -of GTiff "test-translate.tif" 
"test-warp.tif"


Error:
ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do not 
belong to the same celestial body
ERROR 6: Cannot find coordinate operations from 
`PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378.14,298.183263207102,LENGTHUNIT["metre",1,ID["EPSG",9001,PRIMEM["Reference
 
meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122,CONVERSION["Polar
 Stereographic (variant B)",METHOD["Polar Stereographic (variant 
B)",ID["EPSG",9829]],PARAMETER["Latitude of standard 
parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude
 of 
origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False
 easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False 
northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1]]]'
 to `EPSG:4326'


But when the +a and +b in the gdal_translate call are changed, gdalwarp does 
work:
gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 +lon_0=0 
+lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0" -a_ullr 0 -364.11191775 
70.903671186 -4415003.88199764 
"HDF5:\"RAD_NL25_RAC_03H_202401220800.h5\"://image1/image_data" 
"test-translate.tif"

What is the problem?








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


Re: [gdal-dev] Convert alpha to nodata?

2024-01-22 Thread Even Rouault via gdal-dev


Le 22/01/2024 à 14:26, Robert Coup a écrit :

Hi Even,

Thanks for your reply!

On Fri, 19 Jan 2024 at 17:32, Even Rouault 
 wrote:


I can't think an easy way of doing what you want.


Is there a functional gap? /Conceptually/ it feels to me like 
VRTRasterBand.ComplexSource.UseMaskBand should "mark" a pixel to be 
nodata-equivalent (NaN?), and then VRTRasterBand.NoDataValue should 
"write" the nodata-marked pixels with the value.



  ...
  
-1
    
in.tif
1
true
      ...
    
    ...
  


But maybe that's not how it's modelled inside?


The above declaration literally tells "use the mask band of band 1 of 
in.tif and expose it as a Int16 band whose nodatavalue metadata is set 
to -1".  So no "intelligent" things like remapping the values of band 1 
to either their native value or -1 depending on the value of the mask band.





Well, you can definitely do that using a VRT Python pixel
function, that
would essentially do your --calc expression


Thanks for confirming it's possible, I'll have a play.

From what I can see the python pixel functions are only called once,
Not sure what you mean by "called once". Well, they are called each time 
IRasterIO() is called on the VRT band. So if you request 1x1 pixel at a 
time, that's not going to be fast. But if you work on reasonably large 
blocks, the setup cost of going through Python should be amortized.

use numpy and (IIUC) vectorised CPU operations...


Not totally sure what numpy does internally, but there shouldn't be any 
measurable difference between a gdal_calc.py and the equivalent formula 
put as a VRT Python pixel function.



would that make it faster/comparable to one of the C pixel functions 
that's doing a naive array loop over rows & columns?
Potentially... That said if you'd want to use that in a multithreaded 
application to read several VRTs from concurrent threads, then you'd hit 
the Python GIL, so that wouldn't scale well (although it might 
potentially be released by numpy when it goes to native code, so maybe 
that's not so bad). But for single threaded, there's indeed the 
potential that numpy based might be faster than the "naive" C pixel 
functions.



A nice exercice for (Pythonist) contributors would be to add a -f VRT
capability to gdal_calc.py that would essentially automate the
writing
of such VRT using Python pixel functions. Just saying...


*adds to a long list of nice exercises*


Ah ah, this was actually a proposal to the wider audience. I'm sure 
there are lots of people looking for ideas for a first nice contribution 
to GDAL :-)


Even


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


Re: [gdal-dev] Convert alpha to nodata?

2024-01-22 Thread Robert Coup via gdal-dev
Hi Even,

Thanks for your reply!

On Fri, 19 Jan 2024 at 17:32, Even Rouault 
wrote:

> I can't think an easy way of doing what you want.
>

Is there a functional gap? *Conceptually* it feels to me like
VRTRasterBand.ComplexSource.UseMaskBand should "mark" a pixel to be
nodata-equivalent (NaN?), and then VRTRasterBand.NoDataValue should "write"
the nodata-marked pixels with the value.


  ...
  
-1

  in.tif
  1
  true
  ...

...
  


But maybe that's not how it's modelled inside?


> Well, you can definitely do that using a VRT Python pixel function, that
> would essentially do your --calc expression
>

Thanks for confirming it's possible, I'll have a play.

>From what I can see the python pixel functions are only called once,
use numpy and (IIUC) vectorised CPU operations... would that make it
faster/comparable to one of the C pixel functions that's doing a naive
array loop over rows & columns?


> A nice exercice for (Pythonist) contributors would be to add a -f VRT
> capability to gdal_calc.py that would essentially automate the writing
> of such VRT using Python pixel functions. Just saying...


*adds to a long list of nice exercises*

Thanks

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