Issue while creating dataset mask for Landsat 7 ETM+ band data by using rasterio

Question

I am trying to create dataset mask of a particular band of Landsat-7 (B2) & Landsat-8 (TIR) data and then plot the same using the following code:

For Landsat-7 Band-1:-

import rasterio
import matplotlib.pyplot as plt
with rasterio.open(r'D:\Landsat7\B1.tif') as dataset:
    mask = dataset.dataset_mask()
    plt.imshow(mask)

After running the above code I am getting the following plot:

For Landsat-8 Band-1:-

import rasterio
import matplotlib.pyplot as plt
with rasterio.open(r'D:\Landsat8\B1.tif') as dataset:
    mask = dataset.dataset_mask()
    plt.imshow(mask)

After running the above code I am getting the following plot:

Now issue is I want the dataset mask of Landsat-7 to appear in a similar fashion as that of of Landsat-8.

GDAL info of both Landsat - 7 & 8 are as under:

Landsat-8 GDALINFO:

Files: LC08_L1TP_141044_20190324_20190403_01_T1_TIR.tif
Size is 7721, 7871
Coordinate System is:
PROJCS["WGS 84 / UTM zone 45N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32645"]]
Origin = (99885.000000000000000,2676315.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   99885.000, 2676315.000) ( 83d 3'49.31"E, 24d 8'55.20"N)
Lower Left  (   99885.000, 2440185.000) ( 83d 7'31.00"E, 22d 1'14.16"N)
Upper Right (  331515.000, 2676315.000) ( 85d20'28.31"E, 24d11'25.72"N)
Lower Right (  331515.000, 2440185.000) ( 85d22' 1.97"E, 22d 3'29.99"N)
Center      (  215700.000, 2558250.000) ( 84d13'27.74"E, 23d 6'31.09"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
  Overviews: 3861x3936, 1931x1968, 966x984, 483x492, 242x246
  Mask Flags: PER_DATASET
  Overviews of mask band: 3861x3936, 1931x1968, 966x984, 483x492, 242x246

Landsat-7 GDALINFO:

Driver: GTiff/GeoTIFF
Files: LE07_L1TP_147048_20070221_20170105_01_T1_B2.tif
       LE07_L1TP_147048_20070221_20170105_01_T1_B2.aux
       .\LE07_L1TP_147048_20070221_20170105_01_T1_MTL.txt
Size is 7871, 7091
Coordinate System is:
PROJCS["WGS 84 / UTM zone 43N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32643"]]
Origin = (229785.000000000000000,2025315.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
  METADATATYPE=ODL
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  229785.000, 2025315.000) ( 72d26'37.73"E, 18d18' 1.45"N)
Lower Left  (  229785.000, 1812585.000) ( 72d28'13.06"E, 16d22'45.80"N)
Upper Right (  465915.000, 2025315.000) ( 74d40'38.80"E, 18d19' 2.08"N)
Lower Right (  465915.000, 1812585.000) ( 74d40'50.85"E, 16d23'39.70"N)
Center      (  347850.000, 1918950.000) ( 73d34' 4.91"E, 17d21' 3.53"N)
Band 1 Block=7871x1 Type=Byte, ColorInterp=Gray
  Description = Layer_1
  Metadata:
    LAYER_TYPE=athematic

Can someone please help me out in resolving the issue.

@user2856. I have just added gdalinfo for both Landsat - 7 & 8 data. — RRSC NGP, Aug 23 '20 at 03:16

score 2 · Accepted Answer · answered Aug 23 '20 at 12:15

2

Your Landsat 7 dataset has no mask band or nodata value set. You can either create a mask band or set the nodata value.

Assuming nodata = 0

# Without setting nodata
with rasterio.open(r'landsat7.tif') as dataset:
    mask = dataset.dataset_mask()
    print(mask.max(), mask.min())
255 255
Setting nodata (note opening dataset in r+ "edit" mode)
with rasterio.open(r'landsat7.tif', 'r+') as dataset:
    dataset.nodata = 0  # set nodata to 0
    mask = dataset.dataset_mask()
    print(mask.max(), mask.min())
255 0

answered Aug 23 '20 at 12:15

user2856

65,736
6
115
196

Your suggestion worked for me. After incorporating it I was successfully able to plot the dataset mask. – RRSC NGP Aug 23 '20 at 15:52
I have one doubt. After going through your suggestion I am able to successfully plot the dataset mask of L7 raster data. But as you are aware that L7 has scanline error issues and hence the same scanlines are appearing in the dataset mask. Can you suggest how to remove these scanlines that are appearing in the mask, so that it appears same as that of L8 dataset mask. – RRSC NGP Aug 23 '20 at 18:12
gdal_fillnodata perhaps https://gis.stackexchange.com/a/30316/2856 – user2856 Aug 23 '20 at 20:59

Issue while creating dataset mask for Landsat 7 ETM+ band data by using rasterio

1 Answers1

255 255

Setting nodata (note opening dataset in r+ "edit" mode)

255 0