search

inbo / niche_vlaanderen

Python package to run the NICHE Vlaanderen model
https://inbo.github.io/niche_vlaanderen/
MIT License
5 stars 2 forks source link

CRS warning and error due to mini shift in extent? #206

Closed cecileherr closed 4 years ago

cecileherr commented 5 years ago

Niche Python Package question

Description

I get a CRS warning when I try to load the following tifs (here), while it seems to me that they are both in EPSG 31370. When I try to run the model I get an error "error: different grid size or orientation". There might be a minimal shift (< 0.00001 mm) between the 2 (not even visible in ArcGIS, but) QGIS reports the following characteristics:

grid_brasschaat CRS EPSG:31370 - Belge 1972 / Belgian Lambert 72 - Projected Extent 153884.2065000000002328,216245.0000000030559022 : 166975.9605199956858996,227200.0000000000000000 Pixel Size 24.98426339693833143,-24.95444191343267448 (notice they are asymmetric as expected)

bodem CRS EPSG:31370 - Belge 1972 / Belgian Lambert 72 - Projected Extent 153884.2065000000002328,216245.0000000000000000 : 166975.9605199999932665,227200.0000000000000000 Pixel Size 24.9842633969465524,-24.95444191343963425

Could this be the cause of the error? Or am I overlooking something else?

What I Did

import niche_vlaanderen as nv
import matplotlib.pyplot as plt

# A simple model
simple = nv.Niche()

# Set the path for the input files
path = "./input/" 
simple.set_input("mhw", path + "grid_brasschaat.tif")
simple.set_input("mlw", path + "grid_brasschaat.tif")
simple.set_input("soil_code",path + "bodem_v2_ok_new.tif")

# Warning: CRS definitions are not equal!

simple.run(full_model=False, strict_checks=False) 

# error: different grid size or orientation
 ---------------------------------------------------------------------------
SpatialContextError                       Traceback (most recent call last)
<ipython-input-5-03315d4400d0> in <module>()
----> 1 simple.run(full_model=False, strict_checks=False) # an example with a warning for the mean gw levels

~\AppData\Local\Continuum\anaconda3\envs\niche\lib\site-packages\niche_vlaanderen\niche.py in run(self, full_model, deviation, abiotic, strict_checks)
    467                     "Error, different obliged keys are missing")
    468 
--> 469         self._check_input_files(full_model)
    470 
    471         if full_model and not abiotic:

~\AppData\Local\Continuum\anaconda3\envs\niche\lib\site-packages\niche_vlaanderen\niche.py in _check_input_files(self, full_model)
    355             nodata = dst.nodatavals[0]
    356 
--> 357             window = self._context.get_read_window(SpatialContext(dst))
    358             band = dst.read(1, window=window)
    359 

~\AppData\Local\Continuum\anaconda3\envs\niche\lib\site-packages\niche_vlaanderen\spatial_context.py in get_read_window(self, new_sc)
    202         if not self.check_overlap(new_sc):
    203             raise SpatialContextError(
--> 204                 "Error: No overlap between both Spatial contexts."
    205             )
    206 

SpatialContextError: Error: No overlap between both Spatial contexts.
stijnvanhoey commented 5 years ago

When checking both files on their metadata using GDAL: 


> gdalinfo bodem_v2_ok_new.tif
Driver: GTiff/GeoTIFF
Files: bodem_v2_ok_new.tif
          bodem_v2_ok_new.tif.aux.xml
Size is 524, 439
Coordinate System is:
PROJCS["Belge 1972 / Belgian Lambert 72",
    GEOGCS["Belge 1972",
        DATUM["Reseau_National_Belge_1972",
            SPHEROID["International 1924",6378388,297,
                AUTHORITY["EPSG","7022"]],
            TOWGS84[-106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747],
            AUTHORITY["EPSG","6313"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4313"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",51.16666723333333],
    PARAMETER["standard_parallel_2",49.8333339],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",4.367486666666666],
    PARAMETER["false_easting",150000.013],
    PARAMETER["false_northing",5400088.438],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","31370"]]
Origin = (153884.206500000000233,227200.000000000000000)
Pixel Size = (24.984263396946552,-24.954441913439634)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  153884.207,  227200.000) (  4d25'23.67"E, 51d21'18.31"N)
Lower Left  (  153884.207,  216245.000) (  4d25'23.24"E, 51d15'23.85"N)
Upper Right (  166975.961,  227200.000) (  4d36'40.17"E, 51d21'17.46"N)
Lower Right (  166975.961,  216245.000) (  4d36'38.32"E, 51d15'23.00"N)
Center      (  160430.084,  221722.500) (  4d31' 1.35"E, 51d18'20.79"N)
Band 1 Block=524x3 Type=Float32, ColorInterp=Gray
  Min=2.000 Max=15.000 
  Minimum=2.000, Maximum=15.000, Mean=10.842, StdDev=1.505
  NoData Value=-3.4028234663852886e+38
  Metadata:
    STATISTICS_MAXIMUM=15
    STATISTICS_MEAN=10.841776939262
    STATISTICS_MINIMUM=2
    STATISTICS_STDDEV=1.5048828165499
    STATISTICS_VALID_PERCENT=100
> gdalinfo grid_brasschaat.tif
Driver: GTiff/GeoTIFF
Files: grid_brasschaat.tif
       grid_brasschaat.tif.aux.xml
Size is 524, 439
Coordinate System is:
PROJCS["Belge_1972_Belgian_Lambert_72",
    GEOGCS["GCS_Belge 1972",
        DATUM["Reseau_National_Belge_1972",
            SPHEROID["International_1924",6378388,297,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6313"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",51.16666723333333],
    PARAMETER["standard_parallel_2",49.8333339],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",4.367486666666666],
    PARAMETER["false_easting",150000.013],
    PARAMETER["false_northing",5400088.438],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (153884.206500000000233,227200.000000000000000)
Pixel Size = (24.984263396938331,-24.954441913432674)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  153884.207,  227200.000) (  4d25'23.67"E, 51d21'18.31"N)
Lower Left  (  153884.207,  216245.000) (  4d25'23.24"E, 51d15'23.85"N)
Upper Right (  166975.961,  227200.000) (  4d36'40.17"E, 51d21'17.46"N)
Lower Right (  166975.961,  216245.000) (  4d36'38.32"E, 51d15'23.00"N)
Center      (  160430.084,  221722.500) (  4d31' 1.35"E, 51d18'20.79"N)
Band 1 Block=524x1 Type=Float64, ColorInterp=Gray
  Min=0.000 Max=100.000 
  Minimum=0.000, Maximum=100.000, Mean=49.957, StdDev=29.086
  Metadata:
    STATISTICS_MAXIMUM=100
    STATISTICS_MEAN=49.956815454972
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=29.086472440779

There is indeed a small difference in pixel size in between both rasters. I would try to use one of the grids as the reference/mask and match the other values with it using any GIS tooling?

I do not think it would be a good idea to be less strict on the spatial extent matching of the input files.

stijnvanhoey commented 5 years ago

After trying the following adaptation to both rasters using the gdal command line command gdalwarp:

> gdalwarp -tr 25 -25 -tap -t_srs epsg:31370 grid_brasschaat.tif brasschaat.tif
> gdalwarp -tr 25 -25 -tap -t_srs epsg:31370 bodem_v2_ok_new.tif bodem.tif

I can run the following code:

simple = nv.Niche()
simple.set_input("mhw", path + "brasschaat.tif")
simple.set_input("mlw", path + "brasschaat.tif")
simple.set_input("soil_code",path + "bodem.tif")
simple.run(full_model=False, strict_checks=False)

wihout errors on wrong crs or alignments.

cecileherr commented 5 years ago

Ok, many thanks for checking!

I (unfortunately) need to stick to the strange pixel sizes I got from the groundwater modellers, but indeed if I resample everything, it works. I thought there was a tolerance for shifts smaller than 1 cm (or what does this mean exactly?). I do not consider a shift <0.0000001 mm as problematic, but I guess it is not that easy to implement a tolerance for every little source of errors ...

May I suggest to adapt the warning for set_input to "Warning: CRS definitions are not equal or different grid size or orientation are used!" instead of "Warning: CRS definitions are not equal or different "?

johanvdw commented 5 years ago

@cecileherr

It was indeed the intention to ignore shifts <1cm. I remember we tested this before. The comment definitely needs an update.

I will try to figure out why it is not working in this case.

stijnvanhoey commented 5 years ago

@johanvdw make sure to provide the required unit tests as well.

johanvdw commented 5 years ago

@cecileherr The tolerance of 0.01m was not applied everywhere. I made a small adjustment and I think your grid should work now.

You can test by installing the latest version from a conda prompt (before starting python): pip install git+https://github.com/inbo/niche_vlaanderen.git --upgrade [ if git is not installed you will get an error, in that case you should first do conda install git ]

If I'm at INBO I will try to adjust the documentation and add a test. We can then see if we can do an updated release or if we include some other small changes.

Strictly speaking we should perhaps also reconsider the absolute 0.01m rule. If we have a huge grid and a small offset in pixel size this will sum up throughout the grid. It might be better to use 1% of the cell size (cumulative over the grid) as cutoff value. Given the typical study areas we currently have this is more a hypothetical problem.

cecileherr commented 5 years ago

Many thanks for looking at this issue!. Unfortunately I get an error when I try to install the latest version from github:

`Collecting rasterio (from niche-vlaanderen==1.0) Downloading https://files.pythonhosted.org/packages/ce/90/8af39cdfed9a40bd5980 1f101f74eab5f1562d3557806ddfa2ed86d4be0b/rasterio-1.0.22.tar.gz (1.9MB) 100% |¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦¦| 1.9MB 451kB/s Complete output from command python setup.py egg_info: INFO:root:Building on Windows requires extra options to setup.py to locate n eeded GDAL files. More information is available in the README. ERROR: A GDAL API version must be specified. Provide a path to gdal-config u sing a GDAL_CONFIG environment variable or use a GDAL_VERSION environment variab le.

Command "python setup.py egg_info" failed with error code 1 in C:\Users\cecile_h err\AppData\Local\Temp\pip-build-a_bsw4vi\rasterio\ `

johanvdw commented 5 years ago

@cecileherr Ok, my error it seems. You can use: pip install git+https://github.com/inbo/niche_vlaanderen.git (without the upgrade)

cecileherr commented 5 years ago

@johanvdw OK, this new version works with my slightly shifted (< 0.0001 mm) rasters. I get a CRS warning when I load the input rasters (as expected), but I can run the model without errors. Thanks again :-)