Grid size bug?

[hugohartmann]

Describe the bug When importing a WRFnetCDF layer, and visualizing the properties/information of a WRF parameter, I noticed that the Width and Height of the grid is incorrect. I was expecting to see 612x612 for, but the information panel shows me 613x611.

Note: When I open the same file as a raster layer, then the grid size is correct in the information window.

I also wonder if the generated CRS is correct: +proj=lcc +lat_1=50.5 +lat_2=50.5 +lat_0=50.50000381469727 +lon_0=-10 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

Where in the code is this computed? In my namelist, I have the following settings for the lambert projection: truelat1, truelat2 = 50.5, stand_lon=-10, ref_lon=2, ref_lat=50.5. I am curious how this is translated into the Proj.4 representation.

[letmaik]

Grid size is directly taken from the netcdf variable dimensions: https://github.com/GIS4WRF/gis4wrf/blob/686b167d656f5f1835858f1b94de4423091f957c/gis4wrf/core/transforms/wrf_netcdf_to_gdal.py#L134-L137 Have you checked the netcdf file manually with ncdump or similar?

The generated CRS is most likely fine, note that QGIS/GDAL uses a different way of referencing than WRF, therefore the CRS string will typically not have a 1:1 correspondence with all fields in the namelist. The code is here: https://github.com/GIS4WRF/gis4wrf/blob/686b167d656f5f1835858f1b94de4423091f957c/gis4wrf/core/transforms/wrf_netcdf_to_gdal.py#L413-L486

And the functions that generate the proj4 strings are here: https://github.com/GIS4WRF/gis4wrf/blob/686b167d656f5f1835858f1b94de4423091f957c/gis4wrf/core/crs.py#L85-L132

[hugohartmann]

Grid size is directly taken from the netcdf variable dimensions:

gis4wrf/gis4wrf/core/transforms/wrf_netcdf_to_gdal.py

Lines 134 to 137 in 686b167

ds = nc.Dataset(path) try: rows = ds.dimensions['south_north'].size cols = ds.dimensions['west_east'].size

Have you checked the netcdf file manually with ncdump or similar?

Yes, ncdump gives me 612x612. And, as mentioned, when importing the wrfout files as raster file, the grid size is correct as well.

[letmaik]

Can you upload the file here for us to look at it? (Just drag and drop into the comment field if not too huge)

[hugohartmann]

wrfout_3km_20190314_1200_tt_surface.zip

[hugohartmann]

since the grid is rather big, I extracted T2 only, and lat/lon fields required by your code :)

[hugohartmann]

ncdump gives 612x612 on the T2 variable.

[letmaik]

It's a bit mysterious... So, gis4wrf creates a temporary geotiff file when you show a layer, and this file is... tmp.zip and when you run gdalinfo on it, it correctly shows...

Driver: GTiff/GeoTIFF
Size is 612, 612
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",50.5],
    PARAMETER["standard_parallel_2",50.5],
    PARAMETER["latitude_of_origin",50.50000381469727],
    PARAMETER["central_meridian",-10],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-73283.321714020406944,-849469.595121843158267)
Pixel Size = (3007.985985081270428,3001.898629510594219)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -73283.322, -849469.595) ( 10d53'21.95"W, 42d52'34.30"N)
Lower Left  (  -73283.322,  987692.366) ( 11d16'15.96"W, 59d19'47.48"N)
Upper Right ( 1767604.101, -849469.595) ( 10d53' 7.57"E, 40d39' 6.60"N)
Lower Right ( 1767604.101,  987692.366) ( 19d 4'42.34"E, 56d12'52.94"N)
Center      (  847160.390,   69111.386) (  1d59'37.17"E, 50d30'17.63"N)
Band 1 Block=612x3 Type=Float32, ColorInterp=Gray
  Description = 2019-03-14 12:00:00
  NoData Value=32768

If you then open that file (without gis4wrf) in qgis (which uses gdal under the hood), then you can observe the change in dimension sizes. I'm not sure yet if this is just some weird rounding issue and is actually displayed correctly, or if there's a deeper problem. But given that this is related directly to the GeoTIFF backend of QGIS, it probably makes more sense to open an issue with QGIS with the tif file attached.

[hugohartmann]

Thanks, just raised a ticket (https://issues.qgis.org/issues/21581)

[letmaik]

Closing as this is an upstream issue and likely only affects metadata display.