应用开源软件去掉影像中的云

Posted 2021-03-23 gispathfinder

参考《Open Source Geospatial Tools Applications in Earth Observation》—Chapter 17 Case Study on Multispectral Land Cover Classification

Depending on the time of the year and the location on Earth, the satellite image
scenes can contain some clouds. The images used within this case study contain

some cloud cover. Clouds directly affect image quality and can severely hamper the

image analysis and classification. To obtain a scene with minimal clouds, we combine

the downloaded scenes and create a “cloud free” composite. If you are using your own

satellite imagery acquired for another location and it is cloud free than you can skip

the next processing step. The BQA dataset is in a bit-packed format2 and contains

information on water and clouds. This means we need to decode the information

in order to obtain, for instance, a cloud mask. More information on the bit-packing

encoding scheme is explained in Sect. 14.3, where we also implement a tool to create

a cloud mask from the BQA dataset. Alternatively, you can use the official tool, the

L-LDOPE Toolbelt, that is provided by the USGS.3 As an example, we create a cloud

mask with the tool bqa2cloud from Sect. 14.3:

bqa2 cloud -i LC 82070232013160 LGN 00_ BQA .TIF -o LC 82070232013160 LGN 00_ CLD . TIF
bqa 2 cloud -i LC 82070232013192 LGN 00_ BQA . TIF -o LC 82070232013192 LGN 00_ CLD . TIF

We then stack all bands into a single multi-band raster dataset in GeoTIFF format
using the utility gdal_merge.py. The panchromatic band contains little extra
spectral information and the thermal bands have been acquired with a coarser spatial
resolution. Therefore, we do not include bands 8, 10 and 11.

gdal_ merge.py -o LC82070232013160LGN00_ merge.tif -of GTiff -ps
→ 30 30 - separate -co INTERLEAVE = BAND -co COMPRESS = LZW
→ LC82070232013160LGN00_B1.TIF LC82070232013160 LGN00_ B2.TIF
→ LC82070232013160LGN00_B3.TIF LC82070232013160 LGN00_ B4.TIF
→ LC82070232013160LGN00_B5.TIF LC82070232013160 LGN00_ B6.TIF
→ LC82070232013160LGN00_B7.TIF LC82070232013160 LGN00_ B9.TIF
→ LC82070232013160LGN00_CLD.TIF

//

-o LC82070232013160LGN00_merge.tif
      Name of the output multi-band raster dataset.
-of GTiff
      We explicitly set the GeoTIFF image format for the output (not required here as this is the default).
-ps 30 30
     Set a common spatial resolution for all inputs. Most bands are already at this resolution, but the panchromatic band needs to be resampled from 15 to 30 m.
-separate
      Place each input file into a separate stacked band.
-co INTERLEAVE=BAND
     Creation option for band sequential encoding (see Sect. 3.5).
-co COMPRESS=LZW
     Creation option for LZW compressed output (see Sect. 5.5.4).

To create the cloud free composite in Fig. 17.2c based on the two multi-band
raster datasets in Fig. 17.2a, b, we can use the utility pkcomposite from pktools
(see Sect. 12.2).

pkcomposite -i LC82070232013160LGN00. tif -i
→ LC82070232013192LGN00.tif -o LC82070232013_ composite. tif
→ -bndnodata 8 -srcnodata 1 -dstnodata 0 -cr maxndvi -cb 3 -cb 4

//

-i LC82070232013160LGN00.tif -i ...
      Sequence of input raster datasets.
-o LC82070232013_composite.tif
     Name of the output composite based on the two input datasets. myinline-bndnodata 8 Index of the cloud band defining nodata values in the input raster datasets (9th band as the index starts from 0).
-srcnodata 1
     Value indicating nodata, i.e. cloudy values in this case.
-dstnodata 0
     Value written in output raster dataset in case of nodata (i.e. if both scenes are cloudy).
-cr maxndvi
     In case more than one input is cloud free, we select the pixel with the maximum NDVI value.
-cb 3 -cb 4
    The composite bands: the respective indices for the red and near infrared spectral bands in the multi-band input raster to calculate the NDVI on the fly (see Table 17.1)