0

I am having some serious problems attempting to merge the R,G,B bands from Landsat 8 and then apply pan sharpening to them. Unfortunately, I am not familiar enough with the process to figure out where the problem is occurring. I am hoping that providing some details here, I can get some guidance on the correct approach or what I am doing wrong.

First of all, I am using bands B4, B3, B2 for the R,G,B bands, and B8 for the polychromatic band.

Here is the gdalinfo for each:

B4:

Driver: GTiff/GeoTIFF
Files: C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B4.TIF
   C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B4.TIF.ovr
   C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_MTL.txt
Size is 7691, 7811
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 15N",
    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 15N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-93,
            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["unknown"],
        AREA["World - N hemisphere - 96┬░W to 90┬░W - by country"],
        BBOX[0,-96,84,-90]],
    ID["EPSG",32615]]
Data axis to CRS axis mapping: 1,2
Origin = (482085.000000000000000,4581915.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
  METADATATYPE=ODL
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  482085.000, 4581915.000) ( 93d12'51.41"W, 41d23'18.83"N)
Lower Left  (  482085.000, 4347585.000) ( 93d12'27.73"W, 39d16'38.76"N)
Upper Right (  712815.000, 4581915.000) ( 90d27'20.65"W, 41d21'38.26"N)
Lower Right (  712815.000, 4347585.000) ( 90d32' 1.38"W, 39d15' 5.38"N)
Center      (  597450.000, 4464750.000) ( 91d51'10.34"W, 40d19'39.36"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
  Overviews: 2564x2604, 855x868, 285x290, 95x97

B3:

Driver: GTiff/GeoTIFF
Files: C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B3.TIF
       C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B3.TIF.ovr
       C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_MTL.txt
Size is 7691, 7811
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 15N",
    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 15N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-93,
            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["unknown"],
        AREA["World - N hemisphere - 96┬░W to 90┬░W - by country"],
        BBOX[0,-96,84,-90]],
    ID["EPSG",32615]]
Data axis to CRS axis mapping: 1,2
Origin = (482085.000000000000000,4581915.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
  METADATATYPE=ODL
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  482085.000, 4581915.000) ( 93d12'51.41"W, 41d23'18.83"N)
Lower Left  (  482085.000, 4347585.000) ( 93d12'27.73"W, 39d16'38.76"N)
Upper Right (  712815.000, 4581915.000) ( 90d27'20.65"W, 41d21'38.26"N)
Lower Right (  712815.000, 4347585.000) ( 90d32' 1.38"W, 39d15' 5.38"N)
Center      (  597450.000, 4464750.000) ( 91d51'10.34"W, 40d19'39.36"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
  Overviews: 2564x2604, 855x868, 285x290, 95x97 

B2:

Driver: GTiff/GeoTIFF
Files: C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B2.TIF
       C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B2.TIF.ovr
       C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_MTL.txt
Size is 7691, 7811
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 15N",
    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 15N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-93,
            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["unknown"],
        AREA["World - N hemisphere - 96┬░W to 90┬░W - by country"],
        BBOX[0,-96,84,-90]],
    ID["EPSG",32615]]
Data axis to CRS axis mapping: 1,2
Origin = (482085.000000000000000,4581915.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
  METADATATYPE=ODL
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  482085.000, 4581915.000) ( 93d12'51.41"W, 41d23'18.83"N)
Lower Left  (  482085.000, 4347585.000) ( 93d12'27.73"W, 39d16'38.76"N)
Upper Right (  712815.000, 4581915.000) ( 90d27'20.65"W, 41d21'38.26"N)
Lower Right (  712815.000, 4347585.000) ( 90d32' 1.38"W, 39d15' 5.38"N)
Center      (  597450.000, 4464750.000) ( 91d51'10.34"W, 40d19'39.36"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
  Overviews: 2564x2604, 855x868, 285x290, 95x97

B8:

Driver: GTiff/GeoTIFF
Files: C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B8.TIF
       C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_B8.TIF.ovr
       C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\LC08_L1TP_025032_20180916_20180928_01_T1_MTL.txt
Size is 15381, 15621
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 15N",
    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 15N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-93,
            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["unknown"],
        AREA["World - N hemisphere - 96┬░W to 90┬░W - by country"],
        BBOX[0,-96,84,-90]],
    ID["EPSG",32615]]
Data axis to CRS axis mapping: 1,2
Origin = (482092.500000000000000,4581907.500000000000000)
Pixel Size = (15.000000000000000,-15.000000000000000)
Metadata:
  AREA_OR_POINT=Point
  METADATATYPE=ODL
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  482092.500, 4581907.500) ( 93d12'51.08"W, 41d23'18.59"N)
Lower Left  (  482092.500, 4347592.500) ( 93d12'27.42"W, 39d16'39.01"N)
Upper Right (  712807.500, 4581907.500) ( 90d27'20.98"W, 41d21'38.03"N)
Lower Right (  712807.500, 4347592.500) ( 90d32' 1.68"W, 39d15' 5.63"N)
Center      (  597450.000, 4464750.000) ( 91d51'10.34"W, 40d19'39.36"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
  Overviews: 5127x5207, 1709x1736, 570x579, 190x193

I am creating a VRT with the 3 R,G,B bands and then using a Translate to convert it into a RGB Tiff.

Here is the code for the first step:

outds = gdal.BuildVRT(vrt, [f'{r}', f'{g}', f'{b}'], separate=True)
t_options = gdal.TranslateOptions(gdal.ParseCommandLine(f"-co TILED=YES -co PHOTOMETRIC=RGB"))
gdal.Translate(scaled, outds, options=t_options)

where r, g, and b are the paths for the bands 4,3,2 respectively. The corresponding gdalinfo for this image:

Driver: GTiff/GeoTIFF
Files: C:\tmp\aws\Landsat\LC08_L1TP_025032_20180916_20180928_01_T1\im_rescaled.tif
Size is 7691, 7811
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 15N",
    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 15N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-93,
            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["unknown"],
        AREA["World - N hemisphere - 96┬░W to 90┬░W - by country"],
        BBOX[0,-96,84,-90]],
    ID["EPSG",32615]]
Data axis to CRS axis mapping: 1,2
Origin = (482085.000000000000000,4581915.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  482085.000, 4581915.000) ( 93d12'51.41"W, 41d23'18.83"N)
Lower Left  (  482085.000, 4347585.000) ( 93d12'27.73"W, 39d16'38.76"N)
Upper Right (  712815.000, 4581915.000) ( 90d27'20.65"W, 41d21'38.26"N)
Lower Right (  712815.000, 4347585.000) ( 90d32' 1.38"W, 39d15' 5.38"N)
Center      (  597450.000, 4464750.000) ( 91d51'10.34"W, 40d19'39.36"N)
Band 1 Block=256x256 Type=UInt16, ColorInterp=Red
Band 2 Block=256x256 Type=UInt16, ColorInterp=Green
Band 3 Block=256x256 Type=UInt16, ColorInterp=Blue

Now, I can view this in a regular image viewer just fine; albeit the image is very dark. However, if I attempt to view this in ArcGIS Pro or ArcMap, I only receive a white image.

Attempting to pan sharpen via the following code:

command = f"python {basepath}gdal_pansharpen.py {poly} {scaled} {sharpened} -r bilinear -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE"
os.system(command)

Yields a white image in Pro and ArcMap as well. Why can't I view these in ArcMap or Pro?

I have also tried doing a straight pan sharpen using the poly band and the r, g, b bands, but the output for that results in the same all white image. Could that be b/c I haven't applied a scale yet and the values are too high?

I presume I am doing the steps in the correct order? I believe the reason it is so dark is b/c I haven't added any scaling params to the output yet, but I was thinking I had to do that after performing the sharpen?

Can anyone offer any insight? I am using GDAL 3.0.4 for what it's worth.

  • gdal_pansharpen.py *B8.TIF *B4.TIF *B3.TIF *B2.TIF pansharpened_out.tif -r bilinear -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE works for me. – user2856 Mar 16 '20 at 20:31
  • The wildcard may be what's different? – Nathan Raley Mar 17 '20 at 14:52
  • Still end up with the same results... Looks fine in a regular photo viewer, but is all white in a GIS application. – Nathan Raley Mar 17 '20 at 15:20
  • It just spits out all white, does it need to be translated with a scale? – Nathan Raley Mar 17 '20 at 21:19
  • Ah, that looks beautiful. Why'd adding the nodata make it that much better? – Nathan Raley Mar 18 '20 at 16:13
1

All the zeros in your image are throwing out your image statistics so your image appears washed out when stretched for display by your GIS.

Set your NoData value to 0. Alternatively, in your GIS specify 0 as the NoData value / exclude it from the statistics used to stretch (GIS software dependent).

However gdal_pansharpen.py appears to have a bug that you can't set -nodata when specifying separate files as the multispectral bands (#2328):

gdal_pansharpen.py LC08_L1TP_025032_20180916_20180928_01_T1_B8.TIF LC08_L1TP_025032_20180916_20180928_01_T1_B4.TIF LC08_L1TP_025032_20180916_20180928_01_T1_B3.TIF LC08_L1TP_025032_20180916_20180928_01_T1_B2.TIF pansharpened_out.tif -r bilinear -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE -nodata 0
ERROR 5: GDALDataset::GetRasterBand(2) - Illegal band #

So build a VRT:

gdalbuildvrt -separate LC08.vrt LC08_L1TP_025032_20180916_20180928_01_T1_B4.TIF LC08_L1TP_025032_20180916_20180928_01_T1_B3.TIF LC08_L1TP_025032_20180916_20180928_01_T1_B2.TIF
gdal_pansharpen.py -nodata 0 LC08_L1TP_025032_20180916_20180928_01_T1_B8.TIF LC08.vrt pansharpened_out.tif -r bilinear -co PHOTOMETRIC=RGB -co COMPRESS=DEFLATE

enter image description here

  • Ah, that makes it behave much better now. I should probably scale it to the Byte type correct? Should I perform color scale on the pansharpened image statistics or the individual band statistics prior to sharpening? I was trying to follow the post here gis.stackexchange.com/questions/257836/… as mentioned by @radouxju, but could not figure out where he got the last set of numbers from. Thanks! – Nathan Raley Mar 18 '20 at 23:04
  • 1
    Not necessarily... A scaled/stretched RGB byte output is only useful for final visualization. It's not useful for any further analysis and the stretching is fixed, you won't be able to tweak the stretch on the fly in your GIS. All stretch/statistics are from the pan-sharpened image. – user2856 Mar 19 '20 at 0:13
  • We plan on doing analysis on the raw grayscale images, my goal for this is creating a sharp, photo realistic rgb image for the users to see visually. Do you have any idea how @radouxju managed to get the final scale values in his response to the post I listed? The first one, using the mean +/- 1.96 * stddev made sense, I'm not sure how he got his final statistics to get the third set of numbers though? I really appreciate the assist, I pretty much just dove into this a couple of weeks ago. – Nathan Raley Mar 19 '20 at 15:27
  • He's using mean +- 1.96*stdev of a cloud-free subset. Hard to automate the subset selection though. – user2856 Mar 19 '20 at 21:15
  • @NathanRaley Maybe dynamic (based on current extent ) raster stretch could be used here. See: gis.stackexchange.com/questions/257020/… – Comrade Che Mar 20 '20 at 4:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.