search

OSGeo / grass

GRASS GIS - free and open-source geospatial processing engine
https://grass.osgeo.org
Other
852 stars 309 forks source link

[Bug] r.external.out creates raster files with wrong extent/resolution when used with r.tile #3201

Open griembauer opened 1 year ago

griembauer commented 1 year ago

Describe the bug Using r.external.out to directly write to .tif files in combination with r.tile creates result tile raster files that have the correct pixel values but are stretched to match the entire current region extent. Also the pixel size is stretched accordingly. All tile raster files then have the same extent.

To Reproduce In the NC sample location Set the region to the elevation raster

g.region raster=elevation -p
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228500
south:      215000
west:       630000
east:       645000
nsres:      10
ewres:      10
rows:       1350
cols:       1500
cells:      2025000

Run r.external.out

r.external.out directory=/path/to/test/dir extension=".tif" format="GTiff" option="COMPRESS=LZW,PREDICTOR=3,BIGTIFF=YES"

Run r.tile

r.tile input=elevation output=elev_tile width=750 height=675

While the resulting tiles in GRASS have the expected extent, the written raster files are stretched to the region extent. This can be analyzed in a different GIS or by checking any tile in gdalinfo:

gdalinfo /path/to/test/dir/elev_tile-000-000.tif
[...]
Pixel Size = (20.000000000000000,-20.000000000000000)
[...]
Corner Coordinates:
Upper Left  (  630000.000,  228500.000) ( 78d46'27.41"W, 35d48'34.59"N)
Lower Left  (  630000.000,  215000.000) ( 78d46'28.63"W, 35d41'16.54"N)
Upper Right (  645000.000,  228500.000) ( 78d36'29.89"W, 35d48'33.08"N)
Lower Right (  645000.000,  215000.000) ( 78d36'32.01"W, 35d41'15.03"N)
Center      (  637500.000,  221750.000) ( 78d41'29.49"W, 35d44'54.91"N)

(deactivate GDAL output by r.external.out)

r.external.out -rp

Expected behavior Each resulting .tif raster should have the correct extent and resolution such that they form a mosaic of the original input raster used in r.tile

System description (please complete the following information):

metzm commented 11 months ago

r.tile sets initially a single dummy output region for all tiles in L114. It then reads the input raster row by row and distributes the data over the corresponding tiles. This way each input raster row needs to be read only once. Rasters for output tiles are created in L136, using the dummy region which is also used by GDAL to create the output tif. Once an output tile is complete, the correct region is set in L185 with Rast_put_cellhd(). This information is however not passed to GDAL which would need to modify the already set geotransform information. In order to fix this, Rast_put_cellhd() would need to be enabled to update the GDAL output file information. Not sure if this is feasible without substantial changes to the GRASS raster lib and the GRASS GDAL interface therein.