Re: [gdal-dev] HDF5 and geolocation arrays

[Michael Sumner via gdal-dev]

you need

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

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

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

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



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

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

Re: [gdal-dev] Re-post, user-defined

[Liyuneh Gebre via gdal-dev]

Well the arcpy made transformation before re-projection. Appended the
sample script in the shared folder as txt. Anyway how can I achieve this
ArcPy output re-projection from gdal with the given gdalinfo on the input
file  and by editing the PROJ.4 string (target projection )?
PARAMETER["Latitude of natural origin",20,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",6378137,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",5,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],

On Tue, Apr 30, 2024 at 1:11 AM Even Rouault 
wrote:

> At first sight, I'd say that the ArcPy output is incorrect. The SRS
> definition of its output doesn't match what you want since it has
>
> PARAMETER["Latitude of natural origin",20,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8801]],
> PARAMETER["Longitude of natural origin",6378137,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8802]],
> PARAMETER["False easting",5,
> LENGTHUNIT["metre",1],
> ID["EPSG",8806]],
>
> The "Longitude of natural origin" is completely wrong, and should be 5
> given the PROJ.4 string you mentionned
>
> And the the "False easting" should be 0.
>
> The output of GDAL looks plausible.
>
>
> Le 30/04/2024 à 00:00, Liyuneh Gebre a écrit :
>
> Hi Even,
> thank you for your response,
> gdal version is: GDAL 3.6.4, released 2023/04/17
> I got the defined target projection included in the link. yes, it's
> geocentric and for further exploration the custom project transformation
> file .gtf extension is also included.
>
> https://u.pcloud.link/publink/show?code=kZj5wh0ZtFWXtI9UGjQiD6KXxn5A9hiWDpEX
>
> On Mon, Apr 29, 2024 at 4:23 PM Even Rouault 
> wrote:
>
>> Hi,
>>
>> it is difficult to help you with just the elements you've provided.
>>
>> What could be helpful is:
>>
>> - the output of gdalinfo on the input file
>>
>> - the output of gdalinfo on the output file generate by ArcPy
>>
>> - how you define the target projection in ArcPy
>>
>> - the output of gdalinfo on the output file generate by GDAL
>>
>> - links to all the above files
>>
>> - the GDAL version you use
>>
>>
>> Totally blind guess: try to add "+towgs84=0,0,0 " in your -t_srs PROJ
>> string, to ask PROJ to consider that your target CRS datum has the same
>> center and axis orientation as WGS84 (and thus go to geographic <-->
>> geocentric conversion when changing datums). Otherwise without datum
>> shifts, PROJ will consider that longitude,latitude coordinates of the
>> geographic CRS underlying the source and target CRS is the same.
>>
>> Even
>> Le 29/04/2024 à 14:57, Liyuneh Gebre via gdal-dev a écrit :
>>
>> Dear all,
>> I am requesting again if you kind have some responds or direction, thank
>> you
>> I have difficulties in accurately reprojecting a raster file to a
>> user-defined projection.
>> I used the following commands in gdal while the output is not similar to
>> arcpy results.  kindly offer your guidance to achieve the same result.  The
>> arcpy uses custom transformation and then reprojection.
>> gdalwarp -t_srs '+proj=laea +a=6378137.0 +f=0.0 +pm=0 +x_0=0.0 +y_0=0.0
>> +lon_0=20.0 +lat_0=5.0 +units=m +axis=enu +no_defs'  -dstalpha  -wo
>> SOURCE_EXTRA=1000 -r nearest rain_in.bil rain_projected.bil
>>
>>
>> *Liyuneh Gebre**, **Associate Researcher*
>>
>> *Ethiopian Institute of Agricultural Research,EIAR *
>>
>> *Cell Phone:+251 911858155 *
>>
>> *linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
>>  *
>> *skype:liyenew_1*
>> *Addis Ababa, Ethioipa*
>>
>> ___
>> gdal-dev mailing 
>> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> -- http://www.spatialys.com
>> My software is free, but my time generally not.
>>
>>
>
> --
> *Liyuneh Gebre**, **Associate Researcher*
>
> *Ethiopian Institute of Agricultural Research,EIAR *
>
> *Cell Phone:+251 911858155 *
>
> *linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
>  *
> *skype:liyenew_1*
> *Addis Ababa, Ethioipa*
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
*Liyuneh Gebre**, **Associate Researcher*

*Ethiopian Institute of Agricultural Research,EIAR*

*Cell Phone:+251 911858155*

*linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
*
*skype:liyenew_1*
*Addis Ababa, Ethioipa*
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Re-post, user-defined

[Even Rouault via gdal-dev]

At first sight, I'd say that the ArcPy output is incorrect. The SRS 
definition of its output doesn't match what you want since it has


PARAMETER["Latitude of natural origin",20,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",6378137,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",5,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],

The "Longitude of natural origin" is completely wrong, and should be 5 
given the PROJ.4 string you mentionned


And the the "False easting" should be 0.

The output of GDAL looks plausible.


Le 30/04/2024 à 00:00, Liyuneh Gebre a écrit :

Hi Even,
thank you for your response,
gdal version is: GDAL 3.6.4, released 2023/04/17
I got the defined target projection included in the link. yes, it's 
geocentric and for further exploration the custom project 
transformation file .gtf extension is also included.

https://u.pcloud.link/publink/show?code=kZj5wh0ZtFWXtI9UGjQiD6KXxn5A9hiWDpEX

On Mon, Apr 29, 2024 at 4:23 PM Even Rouault 
 wrote:


Hi,

it is difficult to help you with just the elements you've provided.

What could be helpful is:

- the output of gdalinfo on the input file

- the output of gdalinfo on the output file generate by ArcPy

- how you define the target projection in ArcPy

- the output of gdalinfo on the output file generate by GDAL

- links to all the above files

- the GDAL version you use


Totally blind guess: try to add "+towgs84=0,0,0 " in your -t_srs
PROJ string, to ask PROJ to consider that your target CRS datum
has the same center and axis orientation as WGS84 (and thus go to
geographic <--> geocentric conversion when changing datums).
Otherwise without datum shifts, PROJ will consider that
longitude,latitude coordinates of the geographic CRS underlying
the source and target CRS is the same.

Even

Le 29/04/2024 à 14:57, Liyuneh Gebre via gdal-dev a écrit :

Dear all,
I am requesting again if you kind have some responds or
direction, thank you
I have difficulties in accurately reprojecting a raster file to a
user-defined projection.
I used the following commands in gdal while the output is not
similar to arcpy results.  kindly offer your guidance to achieve
the same result.  The arcpy uses custom transformation and then
reprojection.
gdalwarp -t_srs '+proj=laea +a=6378137.0 +f=0.0 +pm=0 +x_0=0.0
+y_0=0.0 +lon_0=20.0 +lat_0=5.0 +units=m +axis=enu +no_defs' 
-dstalpha -wo SOURCE_EXTRA=1000 -r nearest rain_in.bil
rain_projected.bil


/*Liyuneh Gebre*//, Associate Researcher///
/Ethiopian Institute of Agricultural Research,EIAR
/
/Cell Phone:+251 911858155
/
/linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
/
//skype:liyenew_1//
/Addis Ababa, Ethioipa/

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



--
/*Liyuneh Gebre*//, Associate Researcher///
/Ethiopian Institute of Agricultural Research,EIAR
/
/Cell Phone:+251 911858155
/
/linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
/
//skype:liyenew_1//
/Addis Ababa, Ethioipa/


--
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] Re-post, user-defined

[Liyuneh Gebre via gdal-dev]

Hi Even,
thank you for your response,
gdal version is: GDAL 3.6.4, released 2023/04/17
I got the defined target projection included in the link. yes, it's
geocentric and for further exploration the custom project transformation
file .gtf extension is also included.
https://u.pcloud.link/publink/show?code=kZj5wh0ZtFWXtI9UGjQiD6KXxn5A9hiWDpEX

On Mon, Apr 29, 2024 at 4:23 PM Even Rouault 
wrote:

> Hi,
>
> it is difficult to help you with just the elements you've provided.
>
> What could be helpful is:
>
> - the output of gdalinfo on the input file
>
> - the output of gdalinfo on the output file generate by ArcPy
>
> - how you define the target projection in ArcPy
>
> - the output of gdalinfo on the output file generate by GDAL
>
> - links to all the above files
>
> - the GDAL version you use
>
>
> Totally blind guess: try to add "+towgs84=0,0,0 " in your -t_srs PROJ
> string, to ask PROJ to consider that your target CRS datum has the same
> center and axis orientation as WGS84 (and thus go to geographic <-->
> geocentric conversion when changing datums). Otherwise without datum
> shifts, PROJ will consider that longitude,latitude coordinates of the
> geographic CRS underlying the source and target CRS is the same.
>
> Even
> Le 29/04/2024 à 14:57, Liyuneh Gebre via gdal-dev a écrit :
>
> Dear all,
> I am requesting again if you kind have some responds or direction, thank
> you
> I have difficulties in accurately reprojecting a raster file to a
> user-defined projection.
> I used the following commands in gdal while the output is not similar to
> arcpy results.  kindly offer your guidance to achieve the same result.  The
> arcpy uses custom transformation and then reprojection.
> gdalwarp -t_srs '+proj=laea +a=6378137.0 +f=0.0 +pm=0 +x_0=0.0 +y_0=0.0
> +lon_0=20.0 +lat_0=5.0 +units=m +axis=enu +no_defs'  -dstalpha  -wo
> SOURCE_EXTRA=1000 -r nearest rain_in.bil rain_projected.bil
>
>
> *Liyuneh Gebre**, **Associate Researcher*
>
> *Ethiopian Institute of Agricultural Research,EIAR *
>
> *Cell Phone:+251 911858155 *
>
> *linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
>  *
> *skype:liyenew_1*
> *Addis Ababa, Ethioipa*
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
>

-- 
*Liyuneh Gebre**, **Associate Researcher*

*Ethiopian Institute of Agricultural Research,EIAR*

*Cell Phone:+251 911858155*

*linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
*
*skype:liyenew_1*
*Addis Ababa, Ethioipa*
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Test problem After Merge

[Even Rouault via gdal-dev]

Andrew,

Perhaps you're running updated GDAL Python bindings against a libgdal 
that hasn't been rebuilt recently? OSRIsDerivedProjected() has been 
added recently both in libgdal and the bindings


Even

Le 29/04/2024 à 23:27, Andrew Bell via gdal-dev a écrit :

Hi,

I just merged master into my branch and I'm now getting an error when 
trying to run a test. If someone might know what happened, I'd 
appreciate a tip.


(gdal) [master] $ pytest -v autotest/utilities/test_gdal_viewshed.py
ImportError while loading conftest 
'/Users/abell/gdal.2/build/autotest/conftest.py'.

autotest/conftest.py:9: in 
    from osgeo import gdal, ogr, osr
swig/python/osgeo/gdal.py:3942: in 
    from . import ogr
swig/python/osgeo/ogr.py:608: in 
    from . import osr
swig/python/osgeo/osr.py:13: in 
    from . import _osr
E   ImportError: 
dlopen(/Users/abell/gdal.2/build/swig/python/osgeo/_osr.cpython-312-darwin.so 
, 0x0002): symbol not found in flat 
namespace '_OSRIsDerivedProjected'


--
Andrew Bell
andrew.bell...@gmail.com

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


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

Re: [gdal-dev] HDF5 and geolocation arrays

[Joaquim Manuel Freire Luís via gdal-dev]

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

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


Excuse for this little derail, but how do we do that? I mean, the credentials. 
I tried with both:

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

and

gdalinfo  "/vsicurl/https://login:pas...@n5eil01u.ecs.nsidc.org…

but both errored with

ERROR 3: Cannot read 4029 bytes
gdalinfo failed - unable to open 
'/vsicurl/https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5'.

and

ERROR 3: Cannot read 4029 bytes
gdalinfo failed - unable to open 
'/vsicurl/https://login:pas...@n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5'.


Thanks

Joaquim

From: gdal-dev  On Behalf Of Michael Sumner 
via gdal-dev
Sent: Monday, April 29, 2024 8:11 PM
To: gdal-dev 
Subject: [gdal-dev] HDF5 and geolocation arrays

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

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



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

The geolocation arrays are sds 33 and 32 respectively:

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

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

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

I see that the NETCDF driver actually does:

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

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

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

and the two parameter sets are

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

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

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

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

Cheers, Mike


--

Michael Sumner
Software and Database Engineer
Australian Antarctic Division

[gdal-dev] Test problem After Merge

[Andrew Bell via gdal-dev]

Hi,

I just merged master into my branch and I'm now getting an error when
trying to run a test. If someone might know what happened, I'd appreciate a
tip.

(gdal) [master] $ pytest -v autotest/utilities/test_gdal_viewshed.py
ImportError while loading conftest
'/Users/abell/gdal.2/build/autotest/conftest.py'.
autotest/conftest.py:9: in 
from osgeo import gdal, ogr, osr
swig/python/osgeo/gdal.py:3942: in 
from . import ogr
swig/python/osgeo/ogr.py:608: in 
from . import osr
swig/python/osgeo/osr.py:13: in 
from . import _osr
E   ImportError: dlopen(/Users/abell/gdal.2/build/swig/python/osgeo/_
osr.cpython-312-darwin.so, 0x0002): symbol not found in flat namespace
'_OSRIsDerivedProjected'

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

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

[Michael Sumner via gdal-dev]

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

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


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

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

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

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

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

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

Cheers, Mike




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


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

Re: [gdal-dev] HDF5 and geolocation arrays

[Michael Sumner via gdal-dev]

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

Cheers, Mike

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

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

Re: [gdal-dev] HDF5 and geolocation arrays

[Even Rouault via gdal-dev]

Michael,

actually support for those products was already there at 95% since they 
use the HDF5EOS GRID formalism which we support since 3.7. But that 
particular type of products has an oddity. They have an empty 
GROUP=Dimension within the GROUP=GridStructure in the HDF5EOS metadata, 
which our parser didn't like as we use the information in this group to 
map dimension names to dimension index and size. Fixed/worked around in 
https://github.com/OSGeo/gdal/pull/9807 . Hard to tell if those products 
are in fault given that the HDF5EOS specification at 
https://www.earthdata.nasa.gov/s3fs-public/imported/ESDS-RFC-008-v1.1.pdf 
, appendix A, is extremely minimum regarding the HDF5EOS metadata itself 
(for once, I complain for a spec to be too minimal !), and the example 
they give for grids doesn't make sense to me (the declared dimension 
names and sizes have nothing to do with the ones use in the grid)


I now get:

$ gdalinfo 
HDF5:"AMSR_U2_L3_SeaIce12km_B04_20120702.he5"://HDFEOS/GRIDS/SpPolarGrid12km/Data_Fields/SI_12km_SH_18H_ASC

Driver: HDF5Image/HDF5 Dataset
Files: AMSR_U2_L3_SeaIce12km_B04_20120702.he5
Size is 632, 664
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Unknown datum based upon the custom spheroid",
    DATUM["Not specified (based on custom spheroid)",
    ELLIPSOID["Custom spheroid",6378273,0,
    LENGTHUNIT["metre",1,
    ID["EPSG",9001,
    PRIMEM["Greenwich",0,
    ANGLEUNIT["degree",0.0174532925199433,
    ID["EPSG",9122,
    CONVERSION["Polar Stereographic (variant B)",
    METHOD["Polar Stereographic (variant B)",
    ID["EPSG",9829]],
    PARAMETER["Latitude of standard parallel",-70,
    ANGLEUNIT["degree",0.0174532925199433],
    ID["EPSG",8832]],
    PARAMETER["Longitude of origin",0,
    ANGLEUNIT["degree",0.0174532925199433],
    ID["EPSG",8833]],
    PARAMETER["False easting",0,
    LENGTHUNIT["metre",1],
    ID["EPSG",8806]],
    PARAMETER["False northing",0,
    LENGTHUNIT["metre",1],
    ID["EPSG",8807]]],
    CS[Cartesian,2],
    AXIS["(E)",north,
    MERIDIAN[90,
    ANGLEUNIT["degree",0.0174532925199433,
    ID["EPSG",9122]]],
    ORDER[1],
    LENGTHUNIT["Meter",1]],
    AXIS["(N)",north,
    MERIDIAN[0,
    ANGLEUNIT["degree",0.0174532925199433,
    ID["EPSG",9122]]],
    ORDER[2],
    LENGTHUNIT["Meter",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (-395.000,435.000)
Pixel Size = (12500.000,-12500.000)
Metadata:
  Conventions=CF-1.6
  HDFEOS_INFORMATION_HDFEOSVersion=HDFEOS_5.1.15
  history=This version of the Sea Ice processing code contains updates 
provided by the science team on September 16, 2019. For details on these 
updates, see the release notes provided in the DAP.

  institution=NASA's AMSR Science Investigator-led Processing System (SIPS)
  references=Please cite these data as: Markus, T., J. C. Comiso, and 
W. N. Meier. 2018. AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness 
Temperatures, Sea Ice Concentration, Motion & Snow Depth Polar Grids, 
Version 1. [Indicate subset used]. Boulder, Colorado USA. NASA National 
Snow and Ice Data Center Distributed Active Archive Center. doi: 
https://doi.org/10.5067/RA1MIJOYPK3P.

  source=satellite observation
  title=AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures, 
Sea Ice Concentration, Motion & Snow Depth Polar Grids

Corner Coordinates:
Upper Left  (-395.000, 435.000) ( 42d14'27.21"W, 39d11'27.54"S)
Lower Left  (-395.000,-395.000) (135d 0' 0.00"W, 41d23'59.41"S)
Upper Right ( 395.000, 435.000) ( 42d14'27.21"E, 39d11'27.54"S)
Lower Right ( 395.000,-395.000) (135d 0' 0.00"E, 41d23'59.41"S)
Center  (   0.000,  20.000) (  0d 0' 0.01"E, 88d 8'51.76"S)
Band 1 Block=632x1 Type=Int32, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    add_offset=0
    coordinates=lon lat
    long_name=18.7 GHz horizontal daily average ascending Tbs
    packing_convention=netCDF
    packing_convention_description=unpacked = scale_factor x packed + 
add_offset

    scale_factor=0.1
    standard_name=brightness_temperature
    units=degree_kelvin
    _FillValue=0

Even

Le 29/04/2024 à 21:10, Michael Sumner via gdal-dev a écrit :
This HDF5 (requires earthdata credentials your "Authorization: Bearer 
" in GDAL_HTTP_HEADERS, or equiv) presents without geolocation 
arrays.


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

Driver: HDF5Image/HDF5 Dataset
Files: 
/vsicurl/https://n5eil01u.ecs.nsidc.org/AMSA/AU_SI12.001/2012.07.02/AMSR_U2_L3_SeaIce12km_B04_20120702.he5

Size is 608, 896
Metadata:
  Conventions=CF-1.6
  

[gdal-dev] HDF5 and geolocation arrays

[Michael Sumner via gdal-dev]

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

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



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

The geolocation arrays are sds 33 and 32 respectively:

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

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

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

I see that the NETCDF driver actually does:

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

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

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

and the two parameter sets are

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

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

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

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

Cheers, Mike


--

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

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

[Javier Jimenez Shaw via gdal-dev]

Hi all

Thanks for your answers

Just to get some values, the VRT was enough to get an idea. Thanks Jukka.
The questions was more like "this should be somewhere, and I do not find
it", but apparently "it was not there".
I see that GDALRasterInterpolateAtPoint could be useful (I am surprised it
was was not needed before). I will try to add it at some point.

Best regards,
Javier.

On Wed, 24 Apr 2024 at 15:16, Even Rouault 
wrote:

>
> Le 24/04/2024 à 15:00, Michael Sumner a écrit :
>
> Or a grouping function that returned the cell index for neighbours and
> weighting that are involved in whatever calculation summary is wanted.
>
> That doesn't seem super user friendly, as users would still be left to do
> the raster value extraction and applying the weights, taking into account
> nodata, etc. Not trivial. What is the advantage of this compared to
> returning the interpolated value? The only one I see is to potentially save
> a bit of computation if you need to interpolate values at the same location
> in multiple bands, but the performance gain would probably be marginal (or
> if not, then a variant of the function dealing with multiple bands could be
> offered)
>
>
> Maybe the warper could return this as a starting point rather than doing
> the "task at hand". ?
>
> The warper code has indeed a "FilterFuncType
> GWKGetFilterFunc(GDALResampleAlg eResampleAlg)" method that returns a
> function that returns interpolation weights and int
> GWKGetFilterRadius(GDALResampleAlg eResampleAlg). The code in
> GDALRPCGetDEMHeight() has an interesting logic where it caches a window of
> interest around the first queried pixel so that subsequent queries in the
> neighbouroud can be honoured without going to RasterIO(). This
> substantially improves performance in the RPC case, in particular during
> reverse transformation where you use an iterative method and thus may need
> a lot of DEM extraction to compute a single point.
>
>
>
>
> On Wed, Apr 24, 2024 at 8:51 PM Even Rouault via gdal-dev <
> gdal-dev@lists.osgeo.org> wrote:
>
>> Hi,
>>
>> I guess this discussion, and past similar ones, calls for an enhancement.
>> A new API function, like CPLErr
>> GDALRasterInterpolateAtPoint(GDALRasterBandH, double dfPixel, double
>> dfLocation, GDALRIOResampleAlg eInterpolation, double *pdfValue), that
>> could be used by a new mode of gdallocationinfo. The GDALRPCGetDEMHeight()
>> function in alg/gdal_rpc.cpp is a plausible candidate implementation for
>> bilinear and bicubic (we could potentially restrict to that at the moment).
>>
>> Even
>> Le 24/04/2024 à 10:33, Javier Jimenez Shaw via gdal-dev a écrit :
>>
>> Hi
>>
>> I would like to read in QGIS or GDAL an interpolated value in a DSM
>> (well, actually it is a geoid model, but it is the same behaviour). See
>> that I do not want the pixel value, but the linear interpolation among the
>> neighbour pixels, assuming that the pixel value is in the center of the
>> pixel.
>> For instance, this file
>> https://www.isgeoid.polimi.it/Geoid/Asia/Japan/japan2000_g.html
>>
>> Is there any way to get it (without implementing the interpolation
>> myself)?
>> If I understood correctly gdallocationinfo is not interpolating, just
>> giving the pixel value.
>>
>> Thanks
>>
>> .___ ._ ..._ .. . ._.  .___ .. __ . _. . __..  ...  ._ .__
>>
>> ___
>> gdal-dev mailing 
>> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> -- http://www.spatialys.com
>> My software is free, but my time generally not.
>>
>> ___
>> gdal-dev mailing list
>> gdal-dev@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> -- 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] Inquire mapped column names for conversion to shapefile

[Jan Heckman via gdal-dev]

That's pretty cool, you just saved me an awful lot of time!
Thanks,
Jan

On Mon, Apr 29, 2024 at 4:31 PM Even Rouault 
wrote:

> Hi,
>
> Something like:
>
>
> from osgeo import ogr
> original_field_names = [ "toolongforshapefile1", "toolongforshapefile2"]
> map_to_shp = {}
> tmpfilename = "/vsimem/temp.shp"
> ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(tmpfilename)
> lyr = ds.CreateLayer("temp")
> lyr_defn = lyr.GetLayerDefn()
> for name in original_field_names:
>  lyr.CreateField(ogr.FieldDefn(name))
>  shp_field_name = lyr_defn.GetFieldDefn(lyr_defn.GetFieldCount() -
> 1).GetNameRef()
>  map_to_shp[name] = shp_field_name
> ds = None
> ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource(tmpfilename)
> print(map_to_shp)
>
> Even
>
> Le 29/04/2024 à 16:23, Jan Heckman via gdal-dev a écrit :
> > Hi everyone,
> > I want to make a map of original (let's say postgresql/postgis) column
> > names to those used in a shapefile after conversion; assuming that at
> > least some column names in postgis will be over 10 characters in
> > length, and shortening may produce conflicts. Basically I know how
> > this is done, appending _, but I hope that there is preferably a
> > python method available in Qgis python to get the list of column names
> > used in the shapefile.
> > Any help or pointers are welcome!
> > Thanks,
> > Jan
> >
> > ___
> > 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] Inquire mapped column names for conversion to shapefile

[Even Rouault via gdal-dev]

Hi,

Something like:


from osgeo import ogr
original_field_names = [ "toolongforshapefile1", "toolongforshapefile2"]
map_to_shp = {}
tmpfilename = "/vsimem/temp.shp"
ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(tmpfilename)
lyr = ds.CreateLayer("temp")
lyr_defn = lyr.GetLayerDefn()
for name in original_field_names:
    lyr.CreateField(ogr.FieldDefn(name))
    shp_field_name = lyr_defn.GetFieldDefn(lyr_defn.GetFieldCount() - 
1).GetNameRef()

    map_to_shp[name] = shp_field_name
ds = None
ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource(tmpfilename)
print(map_to_shp)

Even

Le 29/04/2024 à 16:23, Jan Heckman via gdal-dev a écrit :

Hi everyone,
I want to make a map of original (let's say postgresql/postgis) column 
names to those used in a shapefile after conversion; assuming that at 
least some column names in postgis will be over 10 characters in 
length, and shortening may produce conflicts. Basically I know how 
this is done, appending _, but I hope that there is preferably a 
python method available in Qgis python to get the list of column names 
used in the shapefile.

Any help or pointers are welcome!
Thanks,
Jan

___
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] Inquire mapped column names for conversion to shapefile

[Jan Heckman via gdal-dev]

Hi everyone,
I want to make a map of original (let's say postgresql/postgis) column
names to those used in a shapefile after conversion; assuming that at least
some column names in postgis will be over 10 characters in length, and
shortening may produce conflicts. Basically I know how this is done,
appending _, but I hope that there is preferably a python method
available in Qgis python to get the list of column names used in the
shapefile.
Any help or pointers are welcome!
Thanks,
Jan
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Re-post, user-defined

[Even Rouault via gdal-dev]

Hi,

it is difficult to help you with just the elements you've provided.

What could be helpful is:

- the output of gdalinfo on the input file

- the output of gdalinfo on the output file generate by ArcPy

- how you define the target projection in ArcPy

- the output of gdalinfo on the output file generate by GDAL

- links to all the above files

- the GDAL version you use


Totally blind guess: try to add "+towgs84=0,0,0 " in your -t_srs PROJ 
string, to ask PROJ to consider that your target CRS datum has the same 
center and axis orientation as WGS84 (and thus go to geographic <--> 
geocentric conversion when changing datums). Otherwise without datum 
shifts, PROJ will consider that longitude,latitude coordinates of the 
geographic CRS underlying the source and target CRS is the same.


Even

Le 29/04/2024 à 14:57, Liyuneh Gebre via gdal-dev a écrit :

Dear all,
I am requesting again if you kind have some responds or direction, 
thank you
I have difficulties in accurately reprojecting a raster file to a 
user-defined projection.
I used the following commands in gdal while the output is not similar 
to arcpy results.  kindly offer your guidance to achieve the same 
result.  The arcpy uses custom transformation and then reprojection.
gdalwarp -t_srs '+proj=laea +a=6378137.0 +f=0.0 +pm=0 +x_0=0.0 
+y_0=0.0 +lon_0=20.0 +lat_0=5.0 +units=m +axis=enu +no_defs'  
-dstalpha -wo SOURCE_EXTRA=1000 -r nearest rain_in.bil rain_projected.bil



/*Liyuneh Gebre*//, Associate Researcher///
/Ethiopian Institute of Agricultural Research,EIAR
/
/Cell Phone:+251 911858155
/
/linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
/
//skype:liyenew_1//
/Addis Ababa, Ethioipa/

___
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] Re-post, user-defined

[Liyuneh Gebre via gdal-dev]

Dear all,
I am requesting again if you kind have some responds or direction, thank you
I have difficulties in accurately reprojecting a raster file to a
user-defined projection.
I used the following commands in gdal while the output is not similar to
arcpy results.  kindly offer your guidance to achieve the same result.  The
arcpy uses custom transformation and then reprojection.
gdalwarp -t_srs '+proj=laea +a=6378137.0 +f=0.0 +pm=0 +x_0=0.0 +y_0=0.0
+lon_0=20.0 +lat_0=5.0 +units=m +axis=enu +no_defs'  -dstalpha  -wo
SOURCE_EXTRA=1000 -r nearest rain_in.bil rain_projected.bil


*Liyuneh Gebre**, **Associate Researcher*

*Ethiopian Institute of Agricultural Research,EIAR*

*Cell Phone:+251 911858155*

*linkedIn: https://www.linkedin.com/in/liyuneh-gebre-b842b440/
*
*skype:liyenew_1*
*Addis Ababa, Ethioipa*
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Reading DXF block attributes

[Simon Gröchenig via gdal-dev]

Hi all,

I attempt to copy DXF data to GPKG using the following command:

ogr2ogr --config DXF_MERGE_BLOCK_GEOMETRIES FALSE --config DXF_INLINE_BLOCKS 
FALSE -f "GPKG" D:/tmp/632_20-2N.gpkg D:/tmp/632_20-2N.dxf --debug ON

My DXF contains some blocks with several block attributes. According to 
https://gdal.org/drivers/vector/dxf.html, the attributes should be written to 
the BlockAttributes field if DXF_INLINE_BLOCKS is set to false. However, this 
field remains empty ("") in my case.

Can you confirm that ogr2ogr can read DXF block attributes? (according to the 
known issues at Github, writing block attributes is not yet supported)

Any ideas why the field BlockAttributes (also BlockScale, BlockOCSNormal, 
BlockOCSCoords) is set to "" in the GPKG? BlockName and BlockAngle are set 
correctly.

Best regards
Simon


___

Vermessungskanzlei Neumayr

Simon Gröchenig - Geoinformation



Albin Egger-Str. 10

9900 Lienz

Tel: +43 4852 68568 31

Mobil: +43 664 1544804

Email: groeche...@zt-gis.at

Web: http://www.zt-gis.at/
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL 3.9.0beta2 available for testing

[Even Rouault via gdal-dev]

Hi,

I've prepared a beta2 of GDAL 3.9.0 as there have been a few non-trivial 
changes since beta1


The major ones are:

* pyproject.toml: use numpy>=2.0.0rc1 for python >=3.9 (#9751)
* Docker: update ubuntu-small and ubuntu-full to 24.04
* TileDB: migrate away from deprecated APIs, and bump minimum version to 
2.15


Other bugfixes:

* CMake: add missing dependency to generate_gdal_version_h target when 
building {plugin}_core targets

* CMake: fix generated find_dependency
* Miramon: various fixes
* GDALSuggestedWarpOutput2(): make changes of PR #9336 honour 
SRC_METHOD=NO_GEOTRANSFORM. Fixes failure on rasterio 
tests/test_warpedvrt.py::test_transformer_options__width_height
* Update p-ranav/argparse to latest master; ogrinfo and ogr2ogr: update 
to use store_into(GIntBig&)

* GDALOpen(): avoid double '.' in error messages related to plugin drivers
* Add OSRIsDerivedProjected() / OGRSpatialReference::IsDerivedProjected()
* OGR_CT: use PROJJSON internally rather than in WKT:2019 (#9732)
* GDALArgumentParser: fix dealing with positional arguments that are 
negative numeric values

* GDALVectorInfo(): fix 3.9.0beta1 regression regarding passing layer names
* Internal libtiff: updated to latest upstream
* gdalinfo and ogrinfo -json: add newline character at end of JSON 
output (qgis#57266)
* OGRFeature::SetField(): add warnings when detecting some lossy 
conversions (refs #9792)
* Arrow/Parquet: fix writing empty point Z with 
GEOMETRY_ENCODING=GEOARROW_INTERLEAVED, and test spatial filtering of 
that encoding (qgis#57228)
* OGRLayer::WriteArrowBatch(): add tolerance for field type mismatches 
if int32/int64/real; Also add an option IF_FIELD_NOT_PRESERVED=ERROR to 
error out when lossy conversion occurs. (#9792)
* Parquet: fix ResetReading() implementation, when using the 
ParquetDataset API and when there's a single batch
* Parquet: fix opening single file Parquet datasets with the 
ParquetDataset API when using PARQUET:filename.parquet

* OGRCloneArrowArray(): add missing support for 'tss:' Arrow type

Docker images with 3.9.0beta2 are currently cooking.

Source snapshots at:

- https://download.osgeo.org/gdal/3.9.0/gdal-3.9.0beta2.tar.gz
- https://download.osgeo.org/gdal/3.9.0/gdal-3.9.0beta2.tar.xz
- https://download.osgeo.org/gdal/3.9.0/gdal390beta2.zip

Autotest snapshots:

- https://download.osgeo.org/gdal/3.9.0/gdalautotest-3.9.0beta2.tar.gz
- https://download.osgeo.org/gdal/3.9.0/gdalautotest-3.9.0beta2.zip

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] Total Viewshed and GDAL

[Rahkonen Jukka via gdal-dev]

Hi,

Can this tool utilize digital surface model and deal with obstacles like 
buildings and trees?

-Jukka Rahkonen-

-Alkuperäinen viesti-
Lähettäjä: gdal-dev  Puolesta Luis Felipe 
Romero via gdal-dev
Lähetetty: maanantai 29. huhtikuuta 2024 12.16
Vastaanottaja: gdal-dev@lists.osgeo.org
Aihe: Re: [gdal-dev] Total Viewshed and GDAL

Hi,

After several months working on other things, I have returned to an application 
of the total viewshed.

What I have done is use Tamas Szekeres' C++ code as a template and, at least on 
Windows, I have managed to compile and run. With just a GPU it is quite fast: 
on a 6000x6000 raster it does the calculations in less than a minute, although 
it can take two or three hours without a GPU on a PC with 8 cores. In multiple 
viewshed mode, the time is reduced proportionally.

On this page I have uploaded more info and some images of a national parks 
project I am working on. https://github.com/luisfromero/ppnn

It's still a little premature to share it, because I need to do many more 
tests, compile in Linux, etc., but my main doubt right now is that I need to 
shape the part of the tool necessary for calculating solar radiation:

In the so-called "mode 3" the visible horizon from each point is calculated, 
which together with other parameters such as tilt and head allow the 
calculation of solar radiation in the presence of shadows.

My question is the file format for the horizon data. For now, it's just an 
array with dimx·dimy·360 byte (using a one degree precision), and I'm 
considering different formats. For example, save it simply as an output file 
(compressed, preferably), but it could be compacted into a geopackage (as a 
metadata file or in a database, I don't know). They could even be saved in a 
multiband raster along with the elevation, tilt and head data from the  DEM. In 
short, everything necessary to then easily calculate the radiation in a fixed 
period of time.

Thanks for any suggestions!

Luis F Romero
___
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] Total Viewshed and GDAL

[Luis Felipe Romero via gdal-dev]

Hi,

After several months working on other things, I have returned to an application 
of the total viewshed.

What I have done is use Tamas Szekeres' C++ code as a template and, at least on Windows, I have managed to compile and run. With just a GPU it is quite fast: on a 6000x6000 raster it does the calculations in less than a minute, although it can take two or three hours without a GPU on a PC with 8 
cores. In multiple viewshed mode, the time is reduced proportionally.


On this page I have uploaded more info and some images of a national parks 
project I am working on. https://github.com/luisfromero/ppnn

It's still a little premature to share it, because I need to do many more 
tests, compile in Linux, etc., but my main doubt right now is that I need to 
shape the part of the tool necessary for calculating solar radiation:

In the so-called "mode 3" the visible horizon from each point is calculated, 
which together with other parameters such as tilt and head allow the calculation of solar 
radiation in the presence of shadows.

My question is the file format for the horizon data. For now, it's just an array with dimx·dimy·360 byte (using a one degree precision), and I'm considering different formats. For example, save it simply as an output file (compressed, preferably), but it could be compacted into a geopackage (as a 
metadata file or in a database, I don't know). They could even be saved in a multiband raster along with the elevation, tilt and head data from the  DEM. In short, everything necessary to then easily calculate the radiation in a fixed period of time.


Thanks for any suggestions!

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

[gdal-dev] Transparent background in hillshade map

[Elena Ruiz via gdal-dev]

Hello, I am using GDAL version 3.6.2 to make a Hillshade map from an elevation 
map (tif) using the gdaldem command.
My question is because I need the background of the image generated with the 
previous command to be transparent,
I have tried different ways but it usually gives me this error:

ERROR 1: Error: band 1 has no color table

These are the commands I have tried:
gdal_translate -b 1 -b 1 -b 1 -b 1 -ot UInt16 --> Generates image with four 
bands, but all in white
gdal_translate -co PHOTOMETRIC=RGB -co ALPHA=YES--> gives ERROR1
gdal_translate -expand rgba --> gives ERROR1
gdalwarp -dstalpha -->creates an image with a transparent background with two 
bands, but AUTOCAD does not recognize it well and considers the entire image 
transparent and not just the background.

Could you give me any solution?
Thanks and regards.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev