Problema durante la creazione della maschera del set di dati per i dati della banda Landsat 7 ETM + utilizzando rasterio
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
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