fadeoutsoftware / WASDI

Web Advanced Space Developers Interface
http://www.wasdi.net
GNU General Public License v3.0
13 stars 5 forks source link

Geoserver hangs with some images (is the NDVI processor faulty?) #477

Open kr1zz opened 2 years ago

kr1zz commented 2 years ago

When visualizing some GeoTIFFs, we get the following symptoms

  1. memory footprint grows quickly
  2. messages like org.geoserver.platform.ServiceException: org.geotools.referencing.operation.projection.ProjectionException: Latitude 5230°15.2'N is too close to a pole. are shown

Sometimes we can see the final image, sometimes tomcat hangs. Why does it happen? Is it a problem of geoserver? Is the tiff we feed broken?

We noticed it happens often with NDVI images, so maybe the S2 indexes processor is worth investigating: is it a problem with SNAP? Is it a problem with GDAL? Is it a problem with our code? There might be a problem with boundaries / projection, although in the end, when the image gets published, it is localized correctly.

Also, the generated image is really big, we should check if we can compress it.

kr1zz commented 2 years ago
$ gdalinfo S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi.tif
Warning 1: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
Driver: GTiff/GeoTIFF
Files: S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi.tif
Size is 10980, 10980
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 33N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 33N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",15,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 12°E to 18°E - by country"],
        BBOX[0,12,84,18]],
    ID["EPSG",32633]]
Data axis to CRS axis mapping: 1,2
Origin = (499980.000000000000000,5400000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_IMAGEDESCRIPTION=S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:

  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  499980.000, 5400000.000) ( 14d59'59.02"E, 48d45'10.85"N)
Lower Left  (  499980.000, 5290200.000) ( 14d59'59.04"E, 47d45'54.60"N)
Upper Right (  609780.000, 5400000.000) ( 16d29'35.77"E, 48d44'36.01"N)
Lower Right (  609780.000, 5290200.000) ( 16d27'53.23"E, 47d45'20.95"N)
Center      (  554880.000, 5345100.000) ( 15d44'21.77"E, 48d15'24.24"N)
Band 1 Block=10980x10980 Type=Float32, ColorInterp=Gray
Band 2 Block=10980x10980 Type=Float32, ColorInterp=Undefined
kr1zz commented 2 years ago

Then I compressed the image:

$ gdal_translate -co "COMPRESS=DEFLATE" S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi.tif S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_DEFLATE.tif

There's a difference in size:

-rw-r--r-- 1 c.nattero 1049089 920M Jan 14 11:41 S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi.tif
-rw-r--r-- 1 c.nattero 1049089 208M Feb  2 10:58 S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_DEFLATE.tif
kr1zz commented 2 years ago
$ gdalinfo S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_DEFLATE.tif
Driver: GTiff/GeoTIFF
Files: S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_DEFLATE.tif
Size is 10980, 10980
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 33N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 33N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",15,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 12°E to 18°E - by country"],
        BBOX[0,12,84,18]],
    ID["EPSG",32633]]
Data axis to CRS axis mapping: 1,2
Origin = (499980.000000000000000,5400000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_IMAGEDESCRIPTION=S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  499980.000, 5400000.000) ( 14d59'59.02"E, 48d45'10.85"N)
Lower Left  (  499980.000, 5290200.000) ( 14d59'59.04"E, 47d45'54.60"N)
Upper Right (  609780.000, 5400000.000) ( 16d29'35.77"E, 48d44'36.01"N)
Lower Right (  609780.000, 5290200.000) ( 16d27'53.23"E, 47d45'20.95"N)
Center      (  554880.000, 5345100.000) ( 15d44'21.77"E, 48d15'24.24"N)
Band 1 Block=10980x1 Type=Float32, ColorInterp=Gray
Band 2 Block=10980x1 Type=Float32, ColorInterp=Undefined
kr1zz commented 2 years ago

I published the compressed version very quickly and without problems

kr1zz commented 2 years ago

Even without compression, gdal_translate seems to solve the issue:

$ gdal_translate S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi.tif S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_UNCOMPRESSED.tif
Warning 1: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
Input file size is 10980, 10980
0...10...20...30...40...50...60...70...80...90...100 - done.

then

$ gdalinfo S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_UNCOMPRESSED.tif
Driver: GTiff/GeoTIFF
Files: S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi_ndvi_UNCOMPRESSED.tif
Size is 10980, 10980
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 33N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 33N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",15,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 12°E to 18°E - by country"],
        BBOX[0,12,84,18]],
    ID["EPSG",32633]]
Data axis to CRS axis mapping: 1,2
Origin = (499980.000000000000000,5400000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_IMAGEDESCRIPTION=S2B_MSIL2A_20220107T095309_N0301_R079_T33UWP_20220107T120820_ndvi
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  499980.000, 5400000.000) ( 14d59'59.02"E, 48d45'10.85"N)
Lower Left  (  499980.000, 5290200.000) ( 14d59'59.04"E, 47d45'54.60"N)
Upper Right (  609780.000, 5400000.000) ( 16d29'35.77"E, 48d44'36.01"N)
Lower Right (  609780.000, 5290200.000) ( 16d27'53.23"E, 47d45'20.95"N)
Center      (  554880.000, 5345100.000) ( 15d44'21.77"E, 48d15'24.24"N)
Band 1 Block=10980x1 Type=Float32, ColorInterp=Gray
Band 2 Block=10980x1 Type=Float32, ColorInterp=Undefined

Uploaded the file, and published the band without problems

ghost commented 2 years ago

If possible, I recommend to test in TEST now we have a standalone Geoserver service :) I'm not saying the standalone Geoserver will fix the problem: I'm pretty sure it won't But we can:

ghost commented 2 years ago

@kr1zz

If you want to retest, now we have a standalone instance everywhere (test + prod):

To know if yes or not we have an improvement :) Maybe it is safer to start in test first :)