mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
717 stars 196 forks source link

coordinate reading bug #116

Closed NengLu closed 4 years ago

NengLu commented 4 years ago

Hi, thanks for developing this good package, but the grid.bbox is not the same as read from gdalinfo, take the dem.tif form this package for example:

!gdalinfo dem.tif
Driver: GTiff/GeoTIFF
Files: dem.tif
Size is 367, 359
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 = (-97.484999999996106,32.821666666665358)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -97.4850000,  32.8216667) ( 97d29' 6.00"W, 32d49'18.00"N)
Lower Left  ( -97.4850000,  32.5225000) ( 97d29' 6.00"W, 32d31'21.00"N)
Upper Right ( -97.1791667,  32.8216667) ( 97d10'45.00"W, 32d49'18.00"N)
Lower Right ( -97.1791667,  32.5225000) ( 97d10'45.00"W, 32d31'21.00"N)
Center      ( -97.3320833,  32.6720833) ( 97d19'55.50"W, 32d40'19.50"N)
Band 1 Block=16x16 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

The grid.bbox should be ( -97.4850000,32.5225000, -97.1791667,32.8216667). But actually:

from pysheds.grid import Grid
grid = Grid.from_raster(dem.tif, data_name='dem')
grid.bbox
(-97.4850000, 32.521667, -97.178333, 32.821667)
mdbartos commented 4 years ago

Greetings,

Thanks for the info. Looks like it's an off-by-one difference.

In GDAL, the lower-right coordinate of the bbox is the upper-left corner of the lower-right pixel.

In pysheds, the lower-right coordinate of the bbox is the lower-right corner of the lower-right pixel:

@property
def bbox(self):
    shape = self.shape
    xmin, ymax = self.affine * (0,0)
    xmax, ymin = self.affine * (shape[1] + 1, shape[0] + 1)
    _bbox = (xmin, ymin, xmax, ymax)
    return _bbox

The desired behavior depends on your application. For instance, the bbox as defined in GDAL doesn't appear to capture the full area of all the cells.

NengLu commented 4 years ago

Thanks for the tips.