Pixel-is-Area vs Pixel-is-Point #2261

CarlosGrohmann · 2022-03-11T14:29:23Z

CarlosGrohmann
Mar 11, 2022
Collaborator

Hi

Just to confirm, which model does GRASS uses? Pixel-is-Area or Pixel-is-Point?

The definitions are in the GeoTiff standard: http://docs.opengeospatial.org/is/19-008r4/19-008r4.html#_pixelisarea_raster_space

"PixelIsArea" Raster Space

The "PixelIsArea" raster grid space R, which is the default, uses coordinates I and J, with (0,0) denoting the upper-left corner of the image, and increasing I to the right, increasing J down. The first pixel-value fills the square grid cell with the bounds:

top-left = (0,0), bottom-right = (1,1)

and so on; by extension this one-by-one grid cell is also referred to as a pixel. An N by M pixel image covers an area with the mathematically defined bounds (0,0), (N,M).


  (0,0) (1,0)
      x---x---x-> I
      | * | * |
      x---x---x       Standard (PixelIsArea) TIFF Raster space R,
      | (1,1) (2,1)   showing the locations (x) and areas (*) of several pixels.
      v
      J

The location of the area pixel is referenced at the upper left of the pixel (pixelCorner convention).

"PixelIsPoint" Raster Space

The PixelIsPoint raster grid space R uses the same coordinate axis names as used in PixelIsArea Raster space, with increasing I to the right, increasing J down. The first pixel-value however, is realized as a point value located at (0,0). An N by M pixel image consists of points which fill the mathematically defined bounds (0,0), (N-1,M-1).

     (0,0)  (1,0)
      x-------x------> I
      |       |
      |       |      PixelIsPoint TIFF Raster space R,
      x-------x      showing the location (x) of several pixels.
      |     (1,1)
      v
      J

The location of the pixel is referenced at the point pixel (pixelCenter convention).

If a point-pixel image were to be displayed on a display device with pixel cells having the same size as the raster spacing, then the upper-left corner of the displayed image would be located in raster space at (-0.5, -0.5).

Thanks

veroandreo · 2022-03-12T20:34:46Z

veroandreo
Mar 12, 2022
Collaborator

Have a look at https://grasswiki.osgeo.org/wiki/GRASS_raster_semantics, I think the answer is there in Cell Locations section

1 reply

CarlosGrohmann Mar 14, 2022
Collaborator Author

Thanks Vero! It is indeed Pixel-is-Area

prof-pguth · 2022-04-02T19:43:37Z

prof-pguth
Apr 2, 2022

In the slope module (r.slope.aspect) I use a Copernicus DEM as the input. It has the pixel-is-point flag (Geotiff key 1025) correctly sect. GRASS changes the grid to pixel-is-area, but does not set Geotiff key 1025 for the output.

Software using the resulting grid is left to guess at the pixel-is-??.

GDAL does the same thing (and gets identical results for the slope grid), but other software winds up with a half pixel offset for the slope grids.

9 replies

hellik Apr 3, 2022
Collaborator

I think the correct solution is for the software (GRASS/GDAL) to set the 1025 GeoKey in the exported grid. The Geotiff specs say:

The use of this geokey is highly recommended for accurate georeferencing of raster.

The uncertainty of not knowing the correct interpretation leads to way too many half-pixel offsets. I think the shift to pixel-is-area is correctly handled, it is just not documented in the output Geotiff.

The way I read the discussion is that GRASS thinks the pixel-is-area is the only logical choice, but the makers of DEMs like DTED, SRTM, or Copernicus never got that message, and the only way for things to be correct is to document the choice in the 1025 Geokey.

.

Regarding SRTM:

https://grass.osgeo.org/grass80/manuals/r.in.srtm.html

AFAIR r.in.srtm correctly interprets the data.

Regarding DTED, there is maybe a wiki entry to set g.region accordingly.

neteler Apr 3, 2022
Maintainer

This issue is GRASS shifting the grid by half a pixel, but then not noting the shift in the Geotiff header.

I have serious doubts that GRASS GIS is shifting the grid.

Please provide a reproducible example.

prof-pguth Apr 3, 2022

Three extracts:

The GRASS command
Geotiff header for the input DEM
Geotiff header for the shifted output slope map

===============================

Use the following command for GRASS:

call "C:\OSGeo4W\bin\o4w_env.bat"
set GDAL_DATA=C:\OSGeo4W\share\gdal
set USE_PATH_FOR_GDAL_PYTHON=YES
grass78 -c c:\mapdata\temp\temp_dem.tif c:\mapdata\temp\grass1\ --exec r.in.gdal input=c:\mapdata\temp\temp_dem.tif output=mymap |more
grass78 c:\mapdata\temp\grass1\PERMANENT --exec r.slope.aspect elevation=mymap slope=slope format=percent |more
grass78 c:\mapdata\temp\grass1\PERMANENT --exec r.out.gdal input=slope out=c:\mapdata\temp\grass_slope_temp_dem.tif type=Float32 --overwrite --quiet |more
rd /S /Q c:\mapdata\temp\grass1

============================================

The Geotiff header for the input file. Note that the corner coordinates are given in field 33932, and both x and y end in 10's. The file also has field 1025 set to indicate RasterPixelIsPoint, which is correct for Copernicus (and SRTM and ASTER).

C:\mapdata\temp\temp_dem.tif
Little Endian (Intel), entries=18

Tag no/TiffTypeName/LengthIm/Offset/TiffTagName

 256    unsigned 2 byte int       1        1387  ImageWidth  
 257    unsigned 2 byte int       1        1749  ImageLength  
 258    unsigned 2 byte int       1          32  BitsPerSample  32
 259    unsigned 2 byte int       1           1  Compression  
 262    unsigned 2 byte int       1           1  PhotometricInterpretation    BlackIsZero
 273    unsigned 2 byte int       1         230  StripOffsets  230
 277    unsigned 2 byte int       1           1  SamplesPerPixel  
 278    unsigned 4 byte int       1        1749  RowsPerStrip  
 279    unsigned 4 byte int       1     9703452  StripByteCounts  
 284    unsigned 2 byte int       1           1  PlanarConfiguration  
 339    unsigned 2 byte int       1           3  Sample format  
 340    double IEEE 8 byte        1     9703682  SMinSampleValue  
 341    double IEEE 8 byte        1     9703690  SMaxSampleValue

33550 double IEEE 8 byte 3 9703698 Geotiff ModelPixelScaleOffset scalex=30 scaley=30 scalez=0
33922 double IEEE 8 byte 6 9703722 Geotiff ModelTiePointOffset Raster coords: 0 0 Model coords: 202020 3197880
34735 unsigned 2 byte int 20 9703770 Geotiff GeoKeyDirectoryOffset
34737 ASCII char 37 9703810 Geotiff GeoASCIIParams Map from MICRODEM UTM||mdz=0/ meters
42113 ASCII char 7 9703847 Missing data (GDAL) -10000

GEOTiff GeoKey Directory Tag
1 Byte 0 4 GeoKeyHeader
1024 other 0 1 1 GTModelTypeGeoKey (6.3.1.1 codes)
1025 other 0 1 2 GTRasterTypeGeoKey (6.3.1.2 codes) RasterPixelIsPoint
1026 other 34737 12 0 GTCitationGeoKey (ASCII)
3072 other 0 1 32628 ProjectedCSTypeGeoKey (6.3.3.1 codes) UTM (ellipse), zone=28N
GEOTIFF registration analyzed OK

=============

Geotiff header for the slope raster created by GRASS. Note that field 33922 has x and y ending in 5's, and that they are in fact shifted by 1/2 pixel (15 m) in each direction. There is no field with the pixel interpretation.

C:\mapdata\temp\grass_slope_temp_dem.tif
Little Endian (Intel), entries=18

Tag no/TiffTypeName/LengthIm/Offset/TiffTagName

 256    unsigned 2 byte int       1        1387  ImageWidth  
 257    unsigned 2 byte int       1        1749  ImageLength  
 258    unsigned 2 byte int       1          32  BitsPerSample  32
 259    unsigned 2 byte int       1           1  Compression  
 262    unsigned 2 byte int       1           1  PhotometricInterpretation    BlackIsZero
 273    unsigned 4 byte int    1749        3728  StripOffsets  8333784,8339332,8344880,8350428,8355976,8361524,8367072,8372620,8378168,8383716,8389264,
 277    unsigned 2 byte int       1           1  SamplesPerPixel  
 278    unsigned 2 byte int       1           1  RowsPerStrip  
 279    unsigned 2 byte int    1749         230  StripByteCounts  5548
 284    unsigned 2 byte int       1           1  PlanarConfiguration  
 305    ASCII char               32       10724  Software  GRASS GIS 7.8.6 with GDAL 3.4.1 
 339    unsigned 2 byte int       1           3  Sample format

33550 double IEEE 8 byte 3 11618 Geotiff ModelPixelScaleOffset scalex=30 scaley=30 scalez=0
33922 double IEEE 8 byte 6 11642 Geotiff ModelTiePointOffset Raster coords: 0 0 Model coords: 202005 3197895
34735 unsigned 2 byte int 32 11690 Geotiff GeoKeyDirectoryOffset
34737 ASCII char 30 11754 Geotiff GeoASCIIParams WGS 84 / UTM zone 28N|WGS 84|
42112 ASCII char 861 10756 Unspecified
7
0.000000e+00 2.000000e+00 255 255 255 255 255 0
2.000000e+00 5.000000e+00 255 255 0 0 255 0
5.000000e+00 1.000000e+01 0 255 0 0 255 255
1.000000e+01 1.500000e+01 0 255 255 0 0 255
1.500000e+01 3.000000e+01 0 0 255 255 0 255
3.000000e+01 5.000000e+01 255 0 255 255 0 0
5.000000e+01 9.000000e+01 255 0 0 0 0 0
slope

42113 ASCII char 4 110 Missing data (GDAL) Õ
GEOTIFF registration analyzed OK

hellik Apr 3, 2022
Collaborator

Three extracts:

The GRASS command

Geotiff header for the input DEM

Geotiff header for the shifted output slope map

===============================

Use the following command for GRASS:

call "C:\OSGeo4W\bin\o4w_env.bat" set GDAL_DATA=C:\OSGeo4W\share\gdal set USE_PATH_FOR_GDAL_PYTHON=YES grass78 -c c:\mapdata\temp\temp_dem.tif c:\mapdata\temp\grass1\ --exec r.in.gdal input=c:\mapdata\temp\temp_dem.tif output=mymap |more grass78 c:\mapdata\temp\grass1\PERMANENT --exec r.slope.aspect elevation=mymap slope=slope format=percent |more grass78 c:\mapdata\temp\grass1\PERMANENT --exec r.out.gdal input=slope out=c:\mapdata\temp\grass_slope_temp_dem.tif type=Float32 --overwrite --quiet |more rd /S /Q c:\mapdata\temp\grass1

============================================

The Geotiff header for the input file. Note that the corner coordinates are given in field 33932, and both x and y end in 10's. The file also has field 1025 set to indicate RasterPixelIsPoint, which is correct for Copernicus (and SRTM and ASTER).

C:\mapdata\temp\temp_dem.tif Little Endian (Intel), entries=18

Tag no/TiffTypeName/LengthIm/Offset/TiffTagName
 256    unsigned 2 byte int       1        1387  ImageWidth  
 257    unsigned 2 byte int       1        1749  ImageLength  
 258    unsigned 2 byte int       1          32  BitsPerSample  32
 259    unsigned 2 byte int       1           1  Compression  
 262    unsigned 2 byte int       1           1  PhotometricInterpretation    BlackIsZero
 273    unsigned 2 byte int       1         230  StripOffsets  230
 277    unsigned 2 byte int       1           1  SamplesPerPixel  
 278    unsigned 4 byte int       1        1749  RowsPerStrip  
 279    unsigned 4 byte int       1     9703452  StripByteCounts  
 284    unsigned 2 byte int       1           1  PlanarConfiguration  
 339    unsigned 2 byte int       1           3  Sample format  
 340    double IEEE 8 byte        1     9703682  SMinSampleValue  
 341    double IEEE 8 byte        1     9703690  SMaxSampleValue  
33550 double IEEE 8 byte 3 9703698 Geotiff ModelPixelScaleOffset scalex=30 scaley=30 scalez=0 33922 double IEEE 8 byte 6 9703722 Geotiff ModelTiePointOffset Raster coords: 0 0 Model coords: 202020 3197880 34735 unsigned 2 byte int 20 9703770 Geotiff GeoKeyDirectoryOffset 34737 ASCII char 37 9703810 Geotiff GeoASCIIParams Map from MICRODEM UTM||mdz=0/ meters 42113 ASCII char 7 9703847 Missing data (GDAL) -10000

GEOTiff GeoKey Directory Tag 1 Byte 0 4 GeoKeyHeader 1024 other 0 1 1 GTModelTypeGeoKey (6.3.1.1 codes) 1025 other 0 1 2 GTRasterTypeGeoKey (6.3.1.2 codes) RasterPixelIsPoint 1026 other 34737 12 0 GTCitationGeoKey (ASCII) 3072 other 0 1 32628 ProjectedCSTypeGeoKey (6.3.3.1 codes) UTM (ellipse), zone=28N GEOTIFF registration analyzed OK

=============

Geotiff header for the slope raster created by GRASS. Note that field 33922 has x and y ending in 5's, and that they are in fact shifted by 1/2 pixel (15 m) in each direction. There is no field with the pixel interpretation.

C:\mapdata\temp\grass_slope_temp_dem.tif Little Endian (Intel), entries=18

Tag no/TiffTypeName/LengthIm/Offset/TiffTagName
 256    unsigned 2 byte int       1        1387  ImageWidth  
 257    unsigned 2 byte int       1        1749  ImageLength  
 258    unsigned 2 byte int       1          32  BitsPerSample  32
 259    unsigned 2 byte int       1           1  Compression  
 262    unsigned 2 byte int       1           1  PhotometricInterpretation    BlackIsZero
 273    unsigned 4 byte int    1749        3728  StripOffsets  8333784,8339332,8344880,8350428,8355976,8361524,8367072,8372620,8378168,8383716,8389264,
 277    unsigned 2 byte int       1           1  SamplesPerPixel  
 278    unsigned 2 byte int       1           1  RowsPerStrip  
 279    unsigned 2 byte int    1749         230  StripByteCounts  5548
 284    unsigned 2 byte int       1           1  PlanarConfiguration  
 305    ASCII char               32       10724  Software  GRASS GIS 7.8.6 with GDAL 3.4.1 
 339    unsigned 2 byte int       1           3  Sample format  
33550 double IEEE 8 byte 3 11618 Geotiff ModelPixelScaleOffset scalex=30 scaley=30 scalez=0 33922 double IEEE 8 byte 6 11642 Geotiff ModelTiePointOffset Raster coords: 0 0 Model coords: 202005 3197895 34735 unsigned 2 byte int 32 11690 Geotiff GeoKeyDirectoryOffset 34737 ASCII char 30 11754 Geotiff GeoASCIIParams WGS 84 / UTM zone 28N|WGS 84| 42112 ASCII char 861 10756 Unspecified 7 0.000000e+00 2.000000e+00 255 255 255 255 255 0 2.000000e+00 5.000000e+00 255 255 0 0 255 0 5.000000e+00 1.000000e+01 0 255 0 0 255 255 1.000000e+01 1.500000e+01 0 255 255 0 0 255 1.500000e+01 3.000000e+01 0 0 255 255 0 255 3.000000e+01 5.000000e+01 255 0 255 255 0 0 5.000000e+01 9.000000e+01 255 0 0 0 0 0 slope

42113 ASCII char 4 110 Missing data (GDAL) Õ GEOTIFF registration analyzed OK

Is this DEM file available?

Which EPSG code has the data?

prof-pguth Apr 3, 2022

Some background:

I am with a group that is trying to compare the global 1 arc second DEMs (SRTM, NASA, ASTER, ALOS, Copernicus). We have published two preliminary papers (https://www.mdpi.com/2072-4292/13/18/3581 and https://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XLIII-B4-2021/395/2021/), the first of which has a discussion in section 4. Point-Based versus Area-Based DEMs.

We want to use slope as a comparison metric, and want to be sure we have reproducible results. The slope profile below highlights that there are three different algorithms in use (Horn, Evans, and Zevenbergen and Thorne) which have relatively minor effect, but that the results from GDAL and GRASS have a horizontal shift compared to MICRODEM and WhiteBox Tools.

We really want to do the slope computations in geographic coordinates so we do not have to deal with the reinterpolation issues (MICRODEM, WhiteBox Tools, and ArcGIS can all do that), but to track down the horizontal offset I started with a UTM DEM which all the programs should be able to handle.

The sample DEM is a reprojected Copernicus DEM, which retains the pixel-is-point interpretation.

The DEM, and slope map from GRASS, are at https://drive.google.com/drive/folders/1n-xEUFAGLBxt9vbDvhLTEW4anvPSAkkC?usp=sharing

The files are EPSG:32628 - WGS 84 / UTM zone 28N

QGIS, MICRODEM, GRASS, GDAL, and WhiteBox Tools all correctly read the projection, but the slope grid does not have an explicit indication that it was changed to Area representation. I have seen references that area is the default, but I cannot find it in the Geotiff specs, which do say it should be explicitly defined.

hellik · 2022-04-03T16:25:56Z

hellik
Apr 3, 2022
Collaborator

System Info                                                                     
GRASS Version: 8.0.1                                                            
Code revision: exported                                                         
Build date: 2022-02-25                                                          
Build platform: x86_64-w64-mingw32                                              
GDAL: 3.4.1                                                                     
PROJ: 8.2.1                                                                     
GEOS: 3.10.2                                                                    
SQLite: 3.37.2                                                                  
Python: 3.9.5                                                                   
wxPython: 4.1.1                                                                 
Platform: Windows-10-10.0.19044-SP0 (OSGeo4W)

gdalinfo temp_dem.tif > dem.txt

Driver: GTiff/GeoTIFF
Files: temp_dem.tif
Size is 1387, 1749
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 28N",
    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 28N",
        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["Engineering survey, topographic mapping."],
        AREA["Between 18°W and 12°W, northern hemisphere between equator and 84°N, onshore and offshore. Gambia. Greenland. Guinea. Guinea-Bissau. Iceland. Ireland - offshore Porcupine Basin. Mauritania. Morocco. Senegal. Sierra Leone. Western Sahara."],
        BBOX[0,-18,84,-12]],
    ID["EPSG",32628]]
Data axis to CRS axis mapping: 1,2
Origin = (202005.000000000000000,3197895.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)
Band 1 Block=1387x1 Type=Float32, ColorInterp=Gray
  NoData Value=-10000

r.in.gdal input=D:\wd\pixel-is-point-issues\temp_dem.tif output=temp_dem        
Importing raster map <temp_dem>...

g.region -p -a raster=temp_dem@data align=temp_dem@data                         
projection: 1 (UTM)
zone:       28
datum:      wgs84
ellipsoid:  wgs84
north:      3197895
south:      3145425
west:       202005
east:       243615
nsres:      30
ewres:      30
rows:       1749
cols:       1387
cells:      2425863

r.info map=temp_dem@data                                                        
 +----------------------------------------------------------------------------+
 | Map:      temp_dem@data                  Date: Sun Apr 03 18:09:42 2022    |
 | Mapset:   data                           Login of Creator: hkmyr           |
 | Location: testdemslope                                                     |
 | DataBase: D:\grassdata                                                     |
 | Title:    temp_dem                                                         |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    FCELL                Semantic label: (none)                |
 |   Rows:         1749                                                       |
 |   Columns:      1387                                                       |
 |   Total Cells:  2425863                                                    |
 |        Projection: UTM (zone 28)                                           |
 |            N:    3197895    S:    3145425   Res:    30                     |
 |            E:     243615    W:     202005   Res:    30                     |
 |   Range of data:    min = -2.515415  max = 2419.683                        |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal input="D:\wd\pixel-is-point-issues\temp_dem.tif" output="t\   |
 |    emp_dem" memory=300 offset=0 num_digits=0                               |
 |                                                                            |
 +----------------------------------------------------------------------------+

r.slope.aspect elevation=temp_dem@data slope=myslope format=percent             
Neigungs-Rasterkarte <myslope> komplett.

r.info map=myslope@data                                                         
 +----------------------------------------------------------------------------+
 | Map:      myslope@data                   Date: Sun Apr 03 18:15:22 2022    |
 | Mapset:   data                           Login of Creator: hkmyr           |
 | Location: testdemslope                                                     |
 | DataBase: D:\grassdata                                                     |
 | Title:    percent slope                                                    |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 385             |
 |   Data Type:    FCELL                Semantic label: (none)                |
 |   Rows:         1749                                                       |
 |   Columns:      1387                                                       |
 |   Total Cells:  2425863                                                    |
 |        Projection: UTM (zone 28)                                           |
 |            N:    3197895    S:    3145425   Res:    30                     |
 |            E:     243615    W:     202005   Res:    30                     |
 |   Range of data:    min = 0  max = 384.0603                                |
 |                                                                            |
 |   Data Source:                                                             |
 |    raster elevation file temp_dem@data                                     |
 |                                                                            |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.slope.aspect                                             |
 |                                                                            |
 |   Comments:                                                                |
 |    slope map elev = temp_dem@data                                          |
 |    zfactor = 1.00 format = percent                                         |
 |    min_slope = 0.000000                                                    |
 |                                                                            |
 |    r.slope.aspect elevation="temp_dem@data" slope="myslope" format="per\   |
 |    cent" precision="FCELL" zscale=1.0 min_slope=0.0                        |
 |                                                                            |
 +----------------------------------------------------------------------------+

r.out.gdal input=myslope@data output=D:\wd\pixel-is-point-issues\myslope.tif format=GTiff
Checking GDAL data type and nodata value...
ERROR 6: D:\wd\pixel-is-point-issues\myslope.tif, band 1:
SetColorTable() only supported for Byte or UInt16 bands in
TIFF format.
Using GDAL data type <Float32>
Die Eingabe-Rasterkarte enthält Zellen mit dem NULL-Wert (no-data). Der Wert nan wird verwendet, um NoData-Werte in der Eingabekarte zu kennzeichnen. Sie können NoData-Wert mit dem Parameter nodata bestimmen.
Exporting raster data to GTiff format...
r.out.gdal komplett. File <D:\wd\pixel-is-point-issues\myslope.tif> created.

D:\wd\pixel-is-point-issues>gdalinfo myslope.tif > myslope.txt

Driver: GTiff/GeoTIFF
Files: myslope.tif
Size is 1387, 1749
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 28N",
    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 28N",
        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["Engineering survey, topographic mapping."],
        AREA["Between 18°W and 12°W, northern hemisphere between equator and 84°N, onshore and offshore. Gambia. Greenland. Guinea. Guinea-Bissau. Iceland. Ireland - offshore Porcupine Basin. Mauritania. Morocco. Senegal. Sierra Leone. Western Sahara."],
        BBOX[0,-18,84,-12]],
    ID["EPSG",32628]]
Data axis to CRS axis mapping: 1,2
Origin = (202005.000000000000000,3197895.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 8.0.1 with GDAL 3.4.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)
Band 1 Block=1387x1 Type=Float32, ColorInterp=Gray
  Description = myslope@data
  NoData Value=nan
  Metadata:
    COLOR_TABLE_RULES_COUNT=7
    COLOR_TABLE_RULE_RGB_0=0.000000e+00 2.000000e+00 255 255 255 255 255 0
    COLOR_TABLE_RULE_RGB_1=2.000000e+00 5.000000e+00 255 255 0 0 255 0
    COLOR_TABLE_RULE_RGB_2=5.000000e+00 1.000000e+01 0 255 0 0 255 255
    COLOR_TABLE_RULE_RGB_3=1.000000e+01 1.500000e+01 0 255 255 0 0 255
    COLOR_TABLE_RULE_RGB_4=1.500000e+01 3.000000e+01 0 0 255 255 0 255
    COLOR_TABLE_RULE_RGB_5=3.000000e+01 5.000000e+01 255 0 255 255 0 0
    COLOR_TABLE_RULE_RGB_6=5.000000e+01 9.000000e+01 255 0 0 0 0 0

0 replies

hellik · 2022-04-03T16:27:59Z

hellik
Apr 3, 2022
Collaborator

extent dem input


Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)

extent slope output

Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 8.0.1 with GDAL 3.4.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)

1 reply

hellik Apr 3, 2022
Collaborator

extent dem input


Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)

extent slope output

Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 8.0.1 with GDAL 3.4.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)

grass_slope_temp_dem.tif
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.8.6 with GDAL 3.4.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)

hellik · 2022-04-03T16:31:42Z

hellik
Apr 3, 2022
Collaborator

D:\wd\pixel-is-point-issues>gdalinfo --version
GDAL 3.4.2, released 2022/03/08

at least gdalinfo reports the same extent of input DEM tif and output slope tif,

0 replies

hellik · 2022-04-03T16:35:35Z

hellik
Apr 3, 2022
Collaborator

and now gdalinfo from the provided slope tif

D:\wd\pixel-is-point-issues>gdalinfo grass_slope_temp_dem.tif
Driver: GTiff/GeoTIFF
Files: grass_slope_temp_dem.tif
Size is 1387, 1749
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 28N",
    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 28N",
        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["Engineering survey, topographic mapping."],
        AREA["Between 18┬░W and 12┬░W, northern hemisphere between equator and 84┬░N, onshore and offshore. Gambia. Greenland. Guinea. Guinea-Bissau. Iceland. Ireland - offshore Porcupine Basin. Mauritania. Morocco. Senegal. Sierra Leone. Western Sahara."],
        BBOX[0,-18,84,-12]],
    ID["EPSG",32628]]
Data axis to CRS axis mapping: 1,2
Origin = (202005.000000000000000,3197895.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.8.6 with GDAL 3.4.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  202005.000, 3197895.000) ( 18d 3'18.11"W, 28d52'27.38"N)
Lower Left  (  202005.000, 3145425.000) ( 18d 2'28.85"W, 28d24' 4.78"N)
Upper Right (  243615.000, 3197895.000) ( 17d37'43.85"W, 28d52'59.76"N)
Lower Right (  243615.000, 3145425.000) ( 17d37' 1.44"W, 28d24'36.53"N)
Center      (  222810.000, 3171660.000) ( 17d50' 8.01"W, 28d38'32.73"N)
Band 1 Block=1387x1 Type=Float32, ColorInterp=Gray
  Description = slope
  NoData Value=nan
  Metadata:
    COLOR_TABLE_RULES_COUNT=7
    COLOR_TABLE_RULE_RGB_0=0.000000e+00 2.000000e+00 255 255 255 255 255 0
    COLOR_TABLE_RULE_RGB_1=2.000000e+00 5.000000e+00 255 255 0 0 255 0
    COLOR_TABLE_RULE_RGB_2=5.000000e+00 1.000000e+01 0 255 0 0 255 255
    COLOR_TABLE_RULE_RGB_3=1.000000e+01 1.500000e+01 0 255 255 0 0 255
    COLOR_TABLE_RULE_RGB_4=1.500000e+01 3.000000e+01 0 0 255 255 0 255
    COLOR_TABLE_RULE_RGB_5=3.000000e+01 5.000000e+01 255 0 255 255 0 0
    COLOR_TABLE_RULE_RGB_6=5.000000e+01 9.000000e+01 255 0 0 0 0 0

0 replies

hellik · 2022-04-03T16:41:11Z

hellik
Apr 3, 2022
Collaborator

Metadata:
  AREA_OR_POINT=Area <=
  TIFFTAG_SOFTWARE=GRASS GIS 8.0.1 with GDAL 3.4.1

0 replies

hellik · 2022-04-03T16:58:03Z

hellik
Apr 3, 2022
Collaborator

according to gdalinfo, the GIS relevant tool, the extend is the same for input and output

7 replies

prof-pguth Apr 3, 2022

I will go there.

It might be relevant for GRASS to point out the change in pixel representation, and its lack of explicit documentation.

hellik Apr 3, 2022
Collaborator

as mentioned above, the output includes this information:

AREA_OR_POINT=Area <=

hellik Apr 3, 2022
Collaborator

it's from the gdalinfo output; so all the necessary GIS information seems to be there.

prof-pguth Apr 3, 2022

Yes, but it should also be in the Geotiff. Forcing anyone who wants to use the Geotiff to use GDAL is counter productive, when GDAL could write a single Geotiff key (and if GDAL only exports Pixel-is-area, very simple coding).

tiffinfo does not appear to be dumping the 4 digit Geotiff tags, which is where this should go.

hellik Apr 3, 2022
Collaborator

That's definitivly a gdal issue

hellik · 2022-04-03T17:11:21Z

hellik
Apr 3, 2022
Collaborator

tiffinfo temp_dem.tif

TIFF Directory at offset 0x8 (8)
  Image Width: 1387 Image Length: 1749
  Bits/Sample: 32
  Sample Format: IEEE floating point
  Compression Scheme: None
  Photometric Interpretation: min-is-black
  Samples/Pixel: 1
  Rows/Strip: 1749
  SMin Sample Value: -2.51541
  SMax Sample Value: 2419.68
  Planar Configuration: single image plane
  Tag 33550: 30.000000,30.000000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,202020.000000,3197880.000000,0.000000
  Tag 34735: 1,1,0,4,1024,0,1,1,1025,0,1,2,1026,34737,12,0,3072,0,1,32628
  Tag 34737: Map from MICRODEM UTM||mdz=0/ meters
  Tag 42113: -10000

tiffinfo grass_slope_temp_dem.tif

TIFF Directory at offset 0x8 (8)
  Image Width: 1387 Image Length: 1749
  Bits/Sample: 32
  Sample Format: IEEE floating point
  Compression Scheme: None
  Photometric Interpretation: min-is-black
  Samples/Pixel: 1
  Rows/Strip: 1
  Planar Configuration: single image plane
  Software: GRASS GIS 7.8.6 with GDAL 3.4.1
  Tag 33550: 30.000000,30.000000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,202005.000000,3197895.000000,0.000000
  Tag 34735: 1,1,0,7,1024,0,1,1,1025,0,1,1,1026,34737,22,0,2049,34737,7,22,2054,0,1,9102,3072,0,1,32628,3076,0,1,9001
  Tag 34737: WGS 84 / UTM zone 28N|WGS 84|
  Tag 42112: <GDALMetadata>
  <Item name="COLOR_TABLE_RULES_COUNT" sample="0">7</Item>
  <Item name="COLOR_TABLE_RULE_RGB_0" sample="0">0.000000e+00 2.000000e+00 255 255 255 255 255 0</Item>
  <Item name="COLOR_TABLE_RULE_RGB_1" sample="0">2.000000e+00 5.000000e+00 255 255 0 0 255 0</Item>
  <Item name="COLOR_TABLE_RULE_RGB_2" sample="0">5.000000e+00 1.000000e+01 0 255 0 0 255 255</Item>
  <Item name="COLOR_TABLE_RULE_RGB_3" sample="0">1.000000e+01 1.500000e+01 0 255 255 0 0 255</Item>
  <Item name="COLOR_TABLE_RULE_RGB_4" sample="0">1.500000e+01 3.000000e+01 0 0 255 255 0 255</Item>
  <Item name="COLOR_TABLE_RULE_RGB_5" sample="0">3.000000e+01 5.000000e+01 255 0 255 255 0 0</Item>
  <Item name="COLOR_TABLE_RULE_RGB_6" sample="0">5.000000e+01 9.000000e+01 255 0 0 0 0 0</Item>
  <Item name="DESCRIPTION" sample="0" role="description">slope</Item>
</GDALMetadata>

  Tag 42113: nan

tiffinfo myslope.tif

TIFF Directory at offset 0x8 (8)
  Image Width: 1387 Image Length: 1749
  Bits/Sample: 32
  Sample Format: IEEE floating point
  Compression Scheme: None
  Photometric Interpretation: min-is-black
  Samples/Pixel: 1
  Rows/Strip: 1
  Planar Configuration: single image plane
  Software: GRASS GIS 8.0.1 with GDAL 3.4.1
  Tag 33550: 30.000000,30.000000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,202005.000000,3197895.000000,0.000000
  Tag 34735: 1,1,0,7,1024,0,1,1,1025,0,1,1,1026,34737,22,0,2049,34737,7,22,2054,0,1,9102,3072,0,1,32628,3076,0,1,9001
  Tag 34737: WGS 84 / UTM zone 28N|WGS 84|
  Tag 42112: <GDALMetadata>
  <Item name="COLOR_TABLE_RULES_COUNT" sample="0">7</Item>
  <Item name="COLOR_TABLE_RULE_RGB_0" sample="0">0.000000e+00 2.000000e+00 255 255 255 255 255 0</Item>
  <Item name="COLOR_TABLE_RULE_RGB_1" sample="0">2.000000e+00 5.000000e+00 255 255 0 0 255 0</Item>
  <Item name="COLOR_TABLE_RULE_RGB_2" sample="0">5.000000e+00 1.000000e+01 0 255 0 0 255 255</Item>
  <Item name="COLOR_TABLE_RULE_RGB_3" sample="0">1.000000e+01 1.500000e+01 0 255 255 0 0 255</Item>
  <Item name="COLOR_TABLE_RULE_RGB_4" sample="0">1.500000e+01 3.000000e+01 0 0 255 255 0 255</Item>
  <Item name="COLOR_TABLE_RULE_RGB_5" sample="0">3.000000e+01 5.000000e+01 255 0 255 255 0 0</Item>
  <Item name="COLOR_TABLE_RULE_RGB_6" sample="0">5.000000e+01 9.000000e+01 255 0 0 0 0 0</Item>
  <Item name="DESCRIPTION" sample="0" role="description">myslope@data</Item>
</GDALMetadata>

  Tag 42113: nan

gdalinfo gives correctly the same GIS extent for input and outpout; AFAIU tiffinfo gives the raw information of the tiffs.

0 replies

prof-pguth · 2022-04-03T17:32:35Z

prof-pguth
Apr 3, 2022

Thanks everyone for looking at this. With some tweaks to my code, and assuming that a Geotiff without the key set will be Pixel-is-area, I can now get the two grids to have identical extent. I have only access to the raw information in the tiffs, and would really like to see GDAL code the pixel representation in the output files.

0 replies

jratike80 · 2022-04-03T20:44:12Z

jratike80
Apr 3, 2022

As far as I know GDAL does set the Pixel-is-area metadata and in GeoTIFF version 1.1 it is mandatory.
Use listgeo and seach for this keyed piece of information:

Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      ...

0 replies

prof-pguth · 2022-04-03T21:48:01Z

prof-pguth
Apr 3, 2022

Thanks. I downloaded a new release of OSGEO4W, which had the same GDAL but a newer GRASS, but it was doing the same thing.

I got listgeo, and found that the Geotiff key was in both the new and old output from GRASS/GDAL. It was not being parsed correctly in the two readers I will using; I have fixed mine, and will report the problem to the author of the other. It's another example of the complexity of Geotiff.

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Pixel-is-Area vs Pixel-is-Point #2261

{{title}}

Replies: 12 comments 18 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Pixel-is-Area vs Pixel-is-Point #2261

CarlosGrohmann Mar 11, 2022 Collaborator

Replies: 12 comments · 18 replies

veroandreo Mar 12, 2022 Collaborator

CarlosGrohmann Mar 14, 2022 Collaborator Author

hellik Apr 3, 2022 Collaborator

neteler Apr 3, 2022 Maintainer

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

hellik Apr 3, 2022 Collaborator

CarlosGrohmann
Mar 11, 2022
Collaborator

Replies: 12 comments 18 replies

veroandreo
Mar 12, 2022
Collaborator

CarlosGrohmann Mar 14, 2022
Collaborator Author

hellik Apr 3, 2022
Collaborator

neteler Apr 3, 2022
Maintainer

hellik Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator

hellik Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator

hellik Apr 3, 2022
Collaborator

hellik Apr 3, 2022
Collaborator

hellik Apr 3, 2022
Collaborator

hellik
Apr 3, 2022
Collaborator