[prev in list] [next in list] [prev in thread] [next in thread] 

List:       grass-user
Subject:    Re: [GRASS-user] raster pixel value
From:       Nikos Alexandris <nik () nikosalexandris ! net>
Date:       2019-03-16 5:59:20
Message-ID: 20190316055920.jp3wy6v6cb5n44oy () imap ! dreamhost ! com
[Download RAW message or body]

* Nikos Alexandris <nik@nikosalexandris.net> [2019-03-14 17:57:16 +0100]:

> * Markus Metz <markus.metz.giswork@gmail.com> [2019-03-14 17:07:49 +0100]:
> 
> > On Thu, Mar 14, 2019 at 3:16 PM Nikos Alexandris <nik@nikosalexandris.net>
> > wrote:
> > > 
> > > Following up, why are there differences between GDAL and GRASS GIS in
> > > the following example?
> > > 
> > > This ftp://ftp.soilgrids.org/data/aggregated/5km/OCDENS_M_sl1_5km_ll.tif
> > > raster map, subject to `gdalinfo`:
> > > ```
> > > gdalinfo OCDENS_M_sl1_5km_ll.tif -nogcp -nomd -norat -noct -nofl
> > > Driver: GTiff/GeoTIFF
> > > Files: OCDENS_M_sl1_5km_ll.tif
> > > Size is 7200, 2987
> > > Coordinate System is:
> > > 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"]]
> > > Origin = (-180.000000000000000,87.370000000000005)
> > > Pixel Size = (0.050000000000000,-0.050000000000000)
> > > Corner Coordinates:
> > > Upper Left  (-180.0000000,  87.3700000) (180d 0' 0.00"W, 87d22'12.00"N)
> > > Lower Left  (-180.0000000, -61.9800000) (180d 0' 0.00"W, 61d58'48.00"S)
> > > Upper Right ( 180.0000000,  87.3700000) (180d 0' 0.00"E, 87d22'12.00"N)
> > > Lower Right ( 180.0000000, -61.9800000) (180d 0' 0.00"E, 61d58'48.00"S)
> > > Center      (   0.0000000,  12.6950000) (  0d 0' 0.00"E, 12d41'42.00"N)
> > > Band 1 Block=7200x1 Type=Int16, ColorInterp=Gray
> > > NoData Value=-32768
> > > ```
> > > 
> > > and GRASS GIS
> > > ```
> > > # create a new Location
> > > grass -c OCDENS_M_sl1_5km_ll.tif /geoyeux/grassdb/global/soil_grids/
> > > 
> > > # projection info
> > > g.proj -g
> > > 
> > > name=WGS 84
> > > datum=wgs84
> > > ellps=wgs84
> > > proj=ll
> > > no_defs=defined
> > > epsg=4326
> > > unit=degree
> > > units=degrees
> > > meters=1.0
> > > 
> > > # raster info
> > > r.info -g OCDENS_M_sl1_5km_ll
> > > 
> > > north=87.37
> > > south=-61.98
> > > east=180
> > > west=-180
> > > nsres=0.05
> > > ewres=0.05
> > > rows=2987
> > > cols=7200
> > > cells=21506400
> > > datatype=CELL
> > > ncats=0
> > > 
> > > # report non-NULL cells and their x, y grid location
> > > r.stats OCDENS_M_sl1_5km_ll -n -x > stats_x
> > > ```
> > > 
> > > Comparing a few single pixels via:
> > > ```
> > > while read LINE;do
> > > set -- $LINE
> > > echo " ($1,$2)
> > > echo "GDAL:  $(gdallocationinfo -valonly OCDENS_M_sl1_5km_ll.tif $1
> > $2)"
> > > echo "GRASS: $3"
> > > echo
> > > done < stats_x_head
> > > ```
> > > 
> > > gives
> > > ```
> > > (2930,77)
> > > GDAL:  2090
> > > GRASS: 2096
> > > 
> > > (2931,77)
> > > GDAL:  2055
> > > GRASS: 2090
> > > 
> > > (2932,77)
> > > GDAL:  2063
> > > GRASS: 2055
> > > 
> > > (2933,77)
> > > GDAL:  2093
> > > GRASS: 2063
> > > 
> > > (2934,77)
> > > GDAL:  2240
> > > GRASS: 2093
> > > 
> > > (2935,77)
> > > GDAL:  2332
> > > GRASS: 2240
> > > 
> > > (2936,77)
> > > GDAL:  2296
> > > GRASS: 2332
> > > 
> > > (2937,77)
> > > GDAL:  2253
> > > GRASS: 2296
> > > 
> > > (2938,77)
> > > GDAL:  2179
> > > GRASS: 2253
> > > 
> > > (2939,77)
> > > GDAL:  2115
> > > GRASS: 2179
> > > ```
> > > 
> > > Why these differences?
> > 
> > With r.stats -x, indexing starts with 1 (first row is 1).
> > With gdallocationinfo, indexing starts with 0 (first row is 0).
> 
> Updated:
> ```
> while read LINE ;do
> set -- ${LINE}
> echo " ($1,$2)"
> echo "GDAL:  $(gdallocationinfo -valonly OCDENS_M_sl1_5km_ll.tif $(echo $1 - 1 |bc) \
> $(echo $2 - 1 |bc))" echo "GRASS: ${3}"
> echo
> done < stats_x_head
> ```

For the sake of clarity, since the file `stats_x_head` is a GRASS GIS
output, where indexing starts from 1:
```
while read LINE ;do
    set -- ${LINE}
    ROW=$(echo $1 -1 |bc)
    COLUMN=$(echo $2 - 1 |bc)
    echo "GDAL  ($ROW, $COLUMN): $(gdallocationinfo -valonly OCDENS_M_sl1_5km_ll.tif \
$ROW $COLUMN)"  echo "GRASS ($1, $2): ${3}"
    echo
done < stats_x_head
```

will return
```

GDAL  (2929, 76): 2096
GRASS (2930, 77): 2096

GDAL  (2930, 76): 2090
GRASS (2931, 77): 2090

GDAL  (2931, 76): 2055
GRASS (2932, 77): 2055

GDAL  (2932, 76): 2063
GRASS (2933, 77): 2063

GDAL  (2933, 76): 2093
GRASS (2934, 77): 2093

GDAL  (2934, 76): 2240
GRASS (2935, 77): 2240

GDAL  (2935, 76): 2332
GRASS (2936, 77): 2332

GDAL  (2936, 76): 2296
GRASS (2937, 77): 2296

GDAL  (2937, 76): 2253
GRASS (2938, 77): 2253

GDAL  (2938, 76): 2179
GRASS (2939, 77): 2179
```

Nikos
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user


[prev in list] [next in list] [prev in thread] [next in thread]