Problema durante la creazione della maschera del set di dati per i dati della banda Landsat 7 ETM + utilizzando rasterio

Aug 22 2020

Sto cercando di creare la maschera del set di dati di una particolare banda di dati Landsat-7 (B2) e Landsat-8 (TIR) ​​e quindi tracciare lo stesso utilizzando il codice seguente:

Per 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)

Dopo aver eseguito il codice sopra, ottengo la seguente trama:

Per 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)

Dopo aver eseguito il codice sopra, ottengo la seguente trama:

Ora il problema è che voglio che la maschera del set di dati del Landsat-7 appaia in modo simile a quella del Landsat-8.

Le informazioni GDAL di entrambi i Landsat - 7 e 8 sono le seguenti:

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

Qualcuno può aiutarmi a risolvere il problema.

Risposte

2 user2856 Aug 23 2020 at 12:15

Il set di dati Landsat 7 non ha una banda di maschera o un valore di nodata impostato. È possibile creare una banda maschera o impostare il valore nodata.

Supponendo 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