[GIS] gdal info hist min and max from python

gdalpython

As pointed at Performing gdal translate to png/jpeg? (by @user30184) gdal info hist gives (STATISTICS_MINIMUM/STATISTICS_MAXIMUN there in the end):

$ gdalinfo -hist image.bsq
Driver: ENVI/ENVI .hdr Labelled
Files: image.bsq
       image.bsq.aux.xml
       image.hdr
Size is 7761, 7891
Coordinate System is:
PROJCS["WGS_1984_UTM_Zone_50N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (318585.000000000000000,4424415.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  Band_1=0
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  318585.000, 4424415.000) (114d52'34.75"E, 39d57' 1.55"N)
Lower Left  (  318585.000, 4187685.000) (114d56'19.99"E, 37d49' 6.80"N)
Upper Right (  551415.000, 4424415.000) (117d36' 7.39"E, 39d58' 5.97"N)
Lower Right (  551415.000, 4187685.000) (117d35' 3.49"E, 37d50' 6.52"N)
Center      (  435000.000, 4306050.000) (116d15' 1.44"E, 38d54' 3.35"N)
Band 1 Block=7761x1 Type=UInt16, ColorInterp=Undefined
  Min=9229.000 Max=20287.000 
  Minimum=9229.000, Maximum=20287.000, Mean=11217.655, StdDev=612.265
  256 buckets from 9207.32 to 20308.7:
  143 372 762 1480 2716 4441 7263 11886 17443 25064 35904 51682 77922 107090 133237 146381 143451 144541 155287 161029 173420 193912 204434 223209 232363 246397 270237 286110 313665 359263 399747 463112 557352 644403 780123 893564 1024421 1174338 1260211 1358161 1473397 1499657 1579140 1569844 1572231 1595593 1535882 1501751 1489654 1381580 1301893 1243696 1114347 1025365 889512 781057 703571 608838 542276 500348 447222 413809 395805 360612 345348 313867 291049 272452 242596 220098 202340 175929 156122 141386 121959 110737 96115 84972 77693 67289 60338 55762 48778 44168 41319 36634 34546 30773 28428 26868 23924 21763 20431 18270 16455 15549 13936 13037 11795 11098 10648 9618 9094 8839 8028 7562 6998 6681 6390 6054 5699 5565 5123 4820 4719 4251 4071 4027 3796 3597 3452 3175 3137 2934 2708 2731 2423 2474 2104 2144 2049 1938 1768 1733 1660 1473 1445 1369 1268 1255 1227 1198 1073 1079 1025 951 872 865 768 792 744 696 711 625 611 584 583 561 512 495 500 459 405 440 403 358 416 335 345 346 347 331 300 301 311 264 276 269 276 258 251 228 217 197 203 221 184 184 176 160 145 149 181 161 157 140 158 168 122 129 115 102 113 101 120 114 98 85 80 86 79 82 68 95 60 67 66 65 55 61 73 78 78 68 43 59 59 71 54 55 60 60 46 41 45 46 33 45 44 40 35 26 36 25 27 20 17 36 24 17 18 17 12 27 14 435 
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=20287
    STATISTICS_MEAN=11217.654590401
    STATISTICS_MINIMUM=9229
    STATISTICS_STDDEV=612.26471631176

I would like to get this values via python so I did:

from osgeo import gdal
info = gdal.Info('image.bsq', deserialize=True)
print(info)

And got:

Driver: ENVI/ENVI .hdr Labelled
Files: image.bsq
       image.bsq.aux.xml
       image.hdr
Size is 7761, 7891
Coordinate System is:
PROJCS["WGS_1984_UTM_Zone_50N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (318585.000000000000000,4424415.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  Band_1=0
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  318585.000, 4424415.000) (114d52'34.75"E, 39d57' 1.55"N)
Lower Left  (  318585.000, 4187685.000) (114d56'19.99"E, 37d49' 6.80"N)
Upper Right (  551415.000, 4424415.000) (117d36' 7.39"E, 39d58' 5.97"N)
Lower Right (  551415.000, 4187685.000) (117d35' 3.49"E, 37d50' 6.52"N)
Center      (  435000.000, 4306050.000) (116d15' 1.44"E, 38d54' 3.35"N)
Band 1 Block=7761x1 Type=UInt16, ColorInterp=Undefined
  Min=9229.000 Max=20287.000 
  Minimum=9229.000, Maximum=20287.000, Mean=11217.655, StdDev=612.265
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=20287
    STATISTICS_MEAN=11217.654590401
    STATISTICS_MINIMUM=9229
    STATISTICS_STDDEV=612.26471631176

How do I get just min/max from all this? So I can:

print(min) to get 9229.

I managet to implement answer by @user30184 but sometimes it does not work. Below is the output of one of this times:

ds = gdal.Open(img_bsq)
band = ds.GetRasterBand(1)
l = band.GetMinimum()
r = band.GetMaximum()
print(gdal.Info(ds))
print(l)
print(r)

From code above:

Driver: ENVI/ENVI .hdr Labelled
Files: B1/image.bsq
       B1/image.bsq.aux.xml
       B1/image.hdr
Size is 7761, 7891
Coordinate System is:
PROJCS["WGS_1984_UTM_Zone_50N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (318585.000000000000000,4424415.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  Band_1=0
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  318585.000, 4424415.000) (114d52'34.75"E, 39d57' 1.55"N)
Lower Left  (  318585.000, 4187685.000) (114d56'19.99"E, 37d49' 6.80"N)
Upper Right (  551415.000, 4424415.000) (117d36' 7.39"E, 39d58' 5.97"N)
Lower Right (  551415.000, 4187685.000) (117d35' 3.49"E, 37d50' 6.52"N)
Center      (  435000.000, 4306050.000) (116d15' 1.44"E, 38d54' 3.35"N)
Band 1 Block=7761x1 Type=UInt16, ColorInterp=Undefined
  NoData Value=0
None
None

Best Answer

Plain GDAL Python solution:

from osgeo import gdal
gtif = gdal.Open( "sentinel.tif" )
srcband = gtif.GetRasterBand(1)
srcband.ComputeStatistics(0)
print "[ MIN ] = ", srcband.GetMinimum()
print "[ MAX ] = ", srcband.GetMaximum()

Comparison of Python and gdalinfo

Python returns:

[ MIN ] =  187.0
[ MAX ] =  18119.0

Gdalinfo returns

  Metadata:
    STATISTICS_MAXIMUM=18119
    STATISTICS_MEAN=3572.9272485136
    STATISTICS_MINIMUM=187
    STATISTICS_STDDEV=1787.6104593137

Related Solutions

[GIS] Rasterize shapefile to a new gtiff file with gdal_rasterize utility,however,I cannot open the tif file

Three suggestions:

Try outputting to another format. Maybe Jpeg or PNG. These will likely work better with Windows Photo Viewer.
Download IrfanView. This will display your tiff quite well.
If neither work - could you run 'gdalinfo -stats Countries02.tif' and post the output.

Edit:

The output is float64. This is likely the issue. Try 'gdal_translate -ot Byte -scale Countries02.tif Countires02_scaled.tif'

I'm not sure why rasterize is creating 64-bit output, but not much reads 64-bit. Also, -scale should not be necessary, but is precautionary.

[GIS] Georeferencing png using gdal

Here's a world file for you (as per Gene's suggestion). Cut and paste the following into a text editor and name the file the same as your full-sized map but give it the extension '.pgw' instead on '.png':

1.0
0.0
0.0
-1.0
0.0
3231.0

I've made the following assumptions:

Your map has a resolution (pixel size) of 1m square (i.e. each pixel is a reasonable stride length for a fantasy game humanoid character). If I am wrong, edit lines 1 and 4 accordingly.
I have assumed that your map is not rotated.
I have assumed that your LL corner is at 0,0 in your world. If your pixel size is bigger than 1m square then you will need to multiply your Y coordinate (3231.0) accordingly.
Technically the UL coordinate (last two lines) is the center of your pixel. I have been lazy and set it to the corner. This in theory means that your coordinates would be consistently half a meter out but it's a fantasy world so nobody will know; it is consistent so won't screw up any distant calculations; and it saves a negative value as one of your origin coordinates. Adjust if you are sufficiently worried :)

I would then re-tile your map from the georeferenced full-sized one (this saves you making a world file for each tile)

Now, if you load your fantasy map into a GIS or similar system, you will probably be asked for the Spatial Reference System (SRS). The World file doesn't encapsulate that. It just locates and orients the image file in space relative to a nominal 0,0 origin. Since this is a fantasy map, all you need to do is pretend that it is in the real world somewhere and pick any coordinate system that gives you a flat projection with units in meters (or feet if you prefer but you will may have to adjust the values for resolution in the world file). A couple of suggestions for SRS could be:

EPSG: 3857 (aka 900913)
EPSG: 27700 (UK National Grid - the advantage of a national grid system UK or otherwise, is that you mostly stay clear of the poles and the international dateline and can hopefully avoid any weirdness which may result from such a location)

Your choice might not matter at all... on the other hand, you choice of SRS might be influenced by what you intend to do with your KML/CSV file (like you might be making a game with a post-apocalyptic theme but set in Quebec ...with zombies...). If you are using Google Earth or similar as a 'platform' for your game/fantasy world, so long as all your other data are in the same coordinate system, all should be well.

However, I do wonder if you may be better sticking with your PNG file as a text format like KML/CSV will be big. A PNG may be better for an internet-based game (reducing bandwidth). If you have a stand-alone scenario, most programming languages understand image formats well (especially if you call the GDAL libraries) and so I'm not convinced that a text-based format is even necessary here - but that's your choice of course.

Best Answer

Related Solutions

[GIS] Rasterize shapefile to a new gtiff file with gdal_rasterize utility,however,I cannot open the tif file

[GIS] Georeferencing png using gdal

Related Question