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.