Re: [Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

2023-09-18 Thread Thomas Lavergne via QGIS-User 
Thank you Even and Richard,

Since Even opened a pull request on gdal's github (
https://github.com/OSGeo/gdal/pull/8407), I will continue the discussion
over there.

Thomas

lør. 16. sep. 2023 kl. 17:10 skrev Even Rouault :

> Thomas,
> >
> > I am not from the GIS world and not even a QGIS user, so bear with me
> > if my answers do not use the expected vocabulary.
> It takes years to get fully familar with GIS oddities, and netCDF is an
> area where even a life long of experience will not be enough to fully
> overcome the "creativity" of the netCDF community...
> >
> > From my point of view, my netCDF files describe their geolocation in a
> > CF-compliant manner, with a grid_mapping CRS variable, and the x and y
> > projection coordinate variables. But QGIS does not recognize the CRS,
> > seemingly because I have my x and y as km (and I do specify
> > :units="km", I am not hiding this).
> >
> > Since I am also a programmer, I tried to find where in the QGIS code
> > the netCDF/CF geolocation is read and decoded, to see if there was a
> > strict test on :units="m". But I failed to find this in the code.
>
> This is actually not really done in QGIS, but in PROJ (more exactly
> QGIS's QgsProjUtils::identifyCrs() delegates to PROJ's proj_identify()
> in
>
> https://github.com/qgis/QGIS/blob/8723b82ec19a3bec7f39b46128c798c7a2ee4230/src/core/proj/qgsprojutils.cpp#L271
> to try to correlate the CRS definition from the netCDF file to one in
> the database).
>
> As you likely open it as a raster (and not a mesh as Richard
> mentionned), the GDAL netCDF driver is also involved since it is it that
> reconstructs a CRS object from the attributes of the netCDF CF conventions.
>
>  From what I can tell, both PROJ and GDAL do the job "as expected". It
> is just that the way this CRS is encoded is not going to find any match
> with a known CRS.  To get what you want, the GDAL netCDF driver should
> both modify the CRS definition to change the km unit to m, and alter the
> geotransform matrix (which gives the coordinate of the top left corner
> and pixel size) to multiply their value by 1000. This could be done, but
> should it be done...? I'm not so sure. The issue is more than the
> "philosophy" of netCDF data producers is sometimes at odds with the
> usual CRS definitions.
>
> That said, for that very particular case, there's a proj4_string
> attribute in Lambert_Azimuthal_Equal_Area variable, with value
> "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m
> +no_defs +type=crs", which indicates the intent to have a metre unit.
> There's also a global attribute geospatial_bounds_crs = "EPSG:6931"
> which further correlate this. Hence this enhancement to the GDAL netCDF
> driver to renormalize to metre: https://github.com/OSGeo/gdal/pull/8407
>
> Even
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
>
>
___
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

Re: [Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

2023-09-16 Thread Even Rouault via QGIS-User 

Thomas,


I am not from the GIS world and not even a QGIS user, so bear with me 
if my answers do not use the expected vocabulary.
It takes years to get fully familar with GIS oddities, and netCDF is an 
area where even a life long of experience will not be enough to fully 
overcome the "creativity" of the netCDF community...


From my point of view, my netCDF files describe their geolocation in a 
CF-compliant manner, with a grid_mapping CRS variable, and the x and y 
projection coordinate variables. But QGIS does not recognize the CRS, 
seemingly because I have my x and y as km (and I do specify 
:units="km", I am not hiding this).


Since I am also a programmer, I tried to find where in the QGIS code 
the netCDF/CF geolocation is read and decoded, to see if there was a 
strict test on :units="m". But I failed to find this in the code.


This is actually not really done in QGIS, but in PROJ (more exactly 
QGIS's QgsProjUtils::identifyCrs() delegates to PROJ's proj_identify() 
in 
https://github.com/qgis/QGIS/blob/8723b82ec19a3bec7f39b46128c798c7a2ee4230/src/core/proj/qgsprojutils.cpp#L271 
to try to correlate the CRS definition from the netCDF file to one in 
the database).


As you likely open it as a raster (and not a mesh as Richard 
mentionned), the GDAL netCDF driver is also involved since it is it that 
reconstructs a CRS object from the attributes of the netCDF CF conventions.


From what I can tell, both PROJ and GDAL do the job "as expected". It 
is just that the way this CRS is encoded is not going to find any match 
with a known CRS.  To get what you want, the GDAL netCDF driver should 
both modify the CRS definition to change the km unit to m, and alter the 
geotransform matrix (which gives the coordinate of the top left corner 
and pixel size) to multiply their value by 1000. This could be done, but 
should it be done...? I'm not so sure. The issue is more than the 
"philosophy" of netCDF data producers is sometimes at odds with the 
usual CRS definitions.


That said, for that very particular case, there's a proj4_string 
attribute in Lambert_Azimuthal_Equal_Area variable, with value 
"+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m 
+no_defs +type=crs", which indicates the intent to have a metre unit. 
There's also a global attribute geospatial_bounds_crs = "EPSG:6931" 
which further correlate this. Hence this enhancement to the GDAL netCDF 
driver to renormalize to metre: https://github.com/OSGeo/gdal/pull/8407


Even

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

___
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

Re: [Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

2023-09-16 Thread Richard Duivenvoorde via QGIS-User 

On 9/16/23 09:37, Thomas Lavergne via QGIS-User wrote:

Since I am also a programmer, I tried to find where in the QGIS code the netCDF/CF 
geolocation is read and decoded, to see if there was a strict test on 
:units="m". But I failed to find this in the code.


Hi Thomas,

I'm a QGIS NetCDF user.

Let me start to say that you can open a netcdf file in QGIS both as a Raster 
and as a Mesh.
And 'netcdf' is a very broad definition with a lot of members, all to be 
handled/defined differently :-)

I think you are using the Mesh dialog, so the Mesh provider:
https://github.com/qgis/QGIS/tree/master/src/providers/mdal

Loading of the crs is done here:
https://github.com/qgis/QGIS/blob/master/src/providers/mdal/qgsmdalprovider.cpp#L493

BUT(!) the actual interpretation (the heavy lifting of all the different Mesh 
types) of the data is done using the MDAL library:
https://github.com/lutraconsulting/MDAL

There are a lot of different mesh data types which MDAL can handle:
https://github.com/lutraconsulting/MDAL/tree/master/mdal/frmts

And as said: all different netcdf flavours have their own handling:
CF files:
https://github.com/lutraconsulting/MDAL/blob/master/mdal/frmts/mdal_cf.hpp
Netcdf ugrids for example:
https://github.com/lutraconsulting/MDAL/blob/master/mdal/frmts/mdal_ugrid.cpp
etc

Not sure if this helps to find out who is to blame:
- the definition/quircks of the data   <= very often :-)
- QGIS
- MDAL

Since you are a programmer, you can try to create your own Netcdf files, for 
example with python
https://unidata.github.io/python-training/workshop/Bonus/netcdf-writing/
and try to create tiny nc files with your exact used dimensions and variables.
Then you can try to find out some scenario's and hopefully find out who is to 
'blame' :-)

(Maybe your actual question is better suited as an issue at the MDAL lib)

Regards,

Richard Duivenvoorde






___
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

Re: [Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

2023-09-16 Thread Thomas Lavergne via QGIS-User 
Thank you Richard,

I am not from the GIS world and not even a QGIS user, so bear with me if my
answers do not use the expected vocabulary.

>From my point of view, my netCDF files describe their geolocation in a
CF-compliant manner, with a grid_mapping CRS variable, and the x and y
projection coordinate variables. But QGIS does not recognize the CRS,
seemingly because I have my x and y as km (and I do specify :units="km", I
am not hiding this).

Since I am also a programmer, I tried to find where in the QGIS code the
netCDF/CF geolocation is read and decoded, to see if there was a strict
test on :units="m". But I failed to find this in the code.

Before I fill out a bug report, I wanted to hear from a QGIS developer (or
someone who knows how to navigate the software code) if we can identify
where the strict test on :units="m" is in the code, and see if alternative
solutions (accepting other units, all multiples of meters) was feasible.

All the best,
Thomas

fre. 15. sep. 2023 kl. 15:16 skrev Richard McDonnell <
richard.mcdonn...@opw.ie>:

> Thomas,
>
> I’m open to correction on this, but meters would generally be the default
> unit used in GIS, when dealing with metric datasets.
>
> What units you use to measure on your Map Canvas or store as attribution
> is then a different thing.
>
>
>
> I had a look for the EPSG codes for EASE2 but all I could get was *EPSG:6931
> <https://epsg.io/6931>*
>
> If you look at the units there you can see meters as the UoM.
>
> I hope that helps,
>
>
>
> Richard.
>
>
>
>
>
>
>
> ——
> Richard McDonnell MSc GIS, FME Certified Professional
> *FRM Data Management*
>
> ——
> Oifig na nOibreacha Poiblí
> Office of Public Works
>
> Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
> Jonathan Swift Street, Trim, Co Meath, C15 NX36
> ——
> M +353 87 688 5964 T +353 46 942 2409
> https://gov.ie/opw
>
> ——
> To send me files larger than 30MB, please use the link below
> https://filetransfer.opw.ie/filedrop/richard.mcdonn...@opw.ie
>
> Email Disclaimer:
> https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/
>
> ——
> MSc GIS, FME Certified Professional
>
> ——
> Oifig na nOibreacha Poiblí
> Office of Public Works
>
> Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
> Jonathan Swift Street, Trim, Co Meath, C15 NX36
> ——
> M +353 87 688 5964 T +353 46 942 2409
> https://https://gov.ie/opw <https://www.opw.ie>
>
> ——
> Email Disclaimer:
> https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/
> <https://www.opw.ie/en/disclaimer/>
>
> *From:* QGIS-User  *On Behalf Of *Thomas
> Lavergne via QGIS-User
> *Sent:* 15 September 2023 13:34
> *To:* qgis-user@lists.osgeo.org
> *Subject:* [Qgis-user] Does QGIS require units "m" for the projection
> coordinate variables in a netCDF/CF file?
>
>
>
> Dear QGIS community,
>
>
>
> We recently released a number of sea-ice climate data records stored in
> netCDF/CF files.
>
> An example file is here
> <https://thredds.met.no/thredds/catalog/osisaf/met.no/reprocessed/ice/drift_455m_files/merged/2020/12/catalog.html?dataset=osisaf/met.no/reprocessed/ice/drift_455m_files/merged/2020/12/ice_drift_nh_ease2-750_cdr-v1p0_24h-202012211200.nc>
> :
>
>
>
> Our files use a polar equal area projection (EASE2) and the x / y axis
> values are given with units "km" (kilometers):
>
>
>
> double xc(xc) ;
> xc:axis = "X" ;
> xc:units = "km" ;
> xc:long_name = "x coordinate of projection (eastings)" ;
> xc:standard_name = "projection_x_coordinate" ;
> double yc(yc) ;
> yc:axis = "Y" ;
> yc:units = "km" ;
> yc:long_name = "y coordinate of projection (northings)" ;
> yc:standard_name = "projection_y_coordinate"
>
>
>
> When opening these in QGIS, they are placed correctly on the map, but the
> message "Unkown CRS" appears in the lower-right corner.
>
>
>
> By manipulating our files, we could get QGIS to recognize the correct CRS.
>
>
>
> ncap2 -s 'xc=xc*1000f;yc=yc*1000f' infile.nc outfile.nc
> ncatted -O -a units,xc,m,c,"m" outfile.nc
> ncatted -O -a units,yc,m,c,"m" outfile.nc
>
>
>
> The NCO commands above change the x / y axis variables to units "m"
> (meters).
>
>
>
> The CF convention
> <http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html>
>  does
> not impose that x / y axis variables are given with unit meters, and shows
> several examples using units "km". But this is apparently an issue for QGIS.
>
>
>
> Can someone please confirm that QGIS requires "meters" for these
> variables, and comment if this is a desired feature or if I should open a
> bug report at github?
>
>
>
> All the best,
>
> Thomas
>
___
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

Re: [Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

2023-09-15 Thread Richard McDonnell via QGIS-User 
Thomas,
I’m open to correction on this, but meters would generally be the default unit 
used in GIS, when dealing with metric datasets.
What units you use to measure on your Map Canvas or store as attribution is 
then a different thing.

I had a look for the EPSG codes for EASE2 but all I could get was 
EPSG:6931<https://epsg.io/6931>
If you look at the units there you can see meters as the UoM.
I hope that helps,

Richard.




——
Richard McDonnell MSc GIS, FME Certified Professional
FRM Data Management

——
Oifig na nOibreacha Poiblí
Office of Public Works

Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
Jonathan Swift Street, Trim, Co Meath, C15 NX36
——
M +353 87 688 5964 T +353 46 942 2409
https://gov.ie/opw

——
To send me files larger than 30MB, please use the link below 
https://filetransfer.opw.ie/filedrop/richard.mcdonn...@opw.ie

Email Disclaimer: 
https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/

——
MSc GIS, FME Certified Professional

——
Oifig na nOibreacha Poiblí
Office of Public Works

Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
Jonathan Swift Street, Trim, Co Meath, C15 NX36
——
M +353 87 688 5964 T +353 46 942 2409
https://https://gov.ie/opw<https://www.opw.ie>

——
Email Disclaimer: 
https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/<https://www.opw.ie/en/disclaimer/>
From: QGIS-User  On Behalf Of Thomas 
Lavergne via QGIS-User
Sent: 15 September 2023 13:34
To: qgis-user@lists.osgeo.org
Subject: [Qgis-user] Does QGIS require units "m" for the projection coordinate 
variables in a netCDF/CF file?

Dear QGIS community,

We recently released a number of sea-ice climate data records stored in 
netCDF/CF files.
An example file is 
here<https://thredds.met.no/thredds/catalog/osisaf/met.no/reprocessed/ice/drift_455m_files/merged/2020/12/catalog.html?dataset=osisaf/met.no/reprocessed/ice/drift_455m_files/merged/2020/12/ice_drift_nh_ease2-750_cdr-v1p0_24h-202012211200.nc>:

Our files use a polar equal area projection (EASE2) and the x / y axis values 
are given with units "km" (kilometers):

double xc(xc) ;
xc:axis = "X" ;
xc:units = "km" ;
xc:long_name = "x coordinate of projection (eastings)" ;
xc:standard_name = "projection_x_coordinate" ;
double yc(yc) ;
yc:axis = "Y" ;
yc:units = "km" ;
yc:long_name = "y coordinate of projection (northings)" ;
yc:standard_name = "projection_y_coordinate"

When opening these in QGIS, they are placed correctly on the map, but the 
message "Unkown CRS" appears in the lower-right corner.

By manipulating our files, we could get QGIS to recognize the correct CRS.

ncap2 -s 'xc=xc*1000f;yc=yc*1000f' infile.nc<http://infile.nc/> 
outfile.nc<http://outfile.nc/>
ncatted -O -a units,xc,m,c,"m" outfile.nc<http://outfile.nc/>
ncatted -O -a units,yc,m,c,"m" outfile.nc<http://outfile.nc/>

The NCO commands above change the x / y axis variables to units "m" (meters).

The CF 
convention<http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html>
 does not impose that x / y axis variables are given with unit meters, and 
shows several examples using units "km". But this is apparently an issue for 
QGIS.

Can someone please confirm that QGIS requires "meters" for these variables, and 
comment if this is a desired feature or if I should open a bug report at github?

All the best,
Thomas
___
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

[Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

2023-09-15 Thread Thomas Lavergne via QGIS-User 
Dear QGIS community,

We recently released a number of sea-ice climate data records stored in
netCDF/CF files.
An example file is here

:

Our files use a polar equal area projection (EASE2) and the x / y axis
values are given with units "km" (kilometers):

double xc(xc) ;
xc:axis = "X" ;
xc:units = "km" ;
xc:long_name = "x coordinate of projection (eastings)" ;
xc:standard_name = "projection_x_coordinate" ;
double yc(yc) ;
yc:axis = "Y" ;
yc:units = "km" ;
yc:long_name = "y coordinate of projection (northings)" ;
yc:standard_name = "projection_y_coordinate"

When opening these in QGIS, they are placed correctly on the map, but the
message "Unkown CRS" appears in the lower-right corner.

By manipulating our files, we could get QGIS to recognize the correct CRS.

ncap2 -s 'xc=xc*1000f;yc=yc*1000f' infile.nc outfile.nc
ncatted -O -a units,xc,m,c,"m" outfile.nc
ncatted -O -a units,yc,m,c,"m" outfile.nc

The NCO commands above change the x / y axis variables to units "m"
(meters).

The CF convention

does
not impose that x / y axis variables are given with unit meters, and shows
several examples using units "km". But this is apparently an issue for QGIS.

Can someone please confirm that QGIS requires "meters" for these variables,
and comment if this is a desired feature or if I should open a bug report
at github?

All the best,
Thomas
___
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user