GIS4WRF / gis4wrf

QGIS toolkit 🧰 for pre- and post-processing 🔨, visualizing 🔍, and running simulations 💻 in the Weather Research and Forecasting (WRF) model 🌀
https://gis4wrf.github.io
MIT License
159 stars 36 forks source link

WPS binary: fix handling of unknown scale, offset, CRS #136

Closed letmaik closed 5 years ago

letmaik commented 5 years ago

Fixes #135.

When scale/offset information is not available for a layer to be converted to WPS binary, then previously this wasn't handled and just raised an assertion error. Now we assume scale=1 offset=0 in that case and carry on. Unfortunately, GDAL's docs are a bit vague on when this may happen ("a boolean to use to indicate if the returned value is meaningful or not"), but it's probably fine...

The other change is when the CRS is missing. In that case we would fail previously, whereas now we simply assume EPSG:4326 (WGS84) as fall-back which will be good enough in those cases. Note that this is a little unfortunate since the user gets prompted by QGIS when opening the layer/file if the CRS is missing, however, the selected CRS is unavailable in our conversion module currently as we work directly on the file level. Forwarding the layer CRS (and other layer modifications like NODATA values) may be a feature to consider. The work-around is that the user saves the layer to a new file (which then has the CRS info) and use that instead. We didn't experience this issue before as we mostly used GeoTIFF files as input which always have full CRS information.

The NetCDF file in #135 exposed all three issues. For reference, this is the ncdump of the original file:

        double FLOATING_POINT_QUANTITY(lat, lon) ;
                ...
                FLOATING_POINT_QUANTITY:coordinates = "lat lon" ;
        float lat(lat) ;
                lat:units = "degrees_north" ;
                lat:standard_name = "latitude" ;
                lat:bounds = "lat_bnds" ;
        float lon(lon) ;
                lon:units = "degrees_east" ;
                lon:standard_name = "longitude" ;
                lon:bounds = "lon_bnds" ;

Re-saving in QGIS after selecting a CRS generated this:

        double Band1(lat, lon) ;
                ...
                Band1:grid_mapping = "crs" ;
        char crs ;
                crs:grid_mapping_name = "latitude_longitude" ;
                crs:long_name = "CRS definition" ;
                crs:longitude_of_prime_meridian = 0. ;
                crs:semi_major_axis = 6378137. ;
                crs:inverse_flattening = 298.257223563 ;
                crs:spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" ;
                crs:GeoTransform = "90 1 0 0 0 -1 " ;
        double lat(lat) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
        double lon(lon) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;