用于批处理 Sentinel-1 图像的 python 脚本

该代码实现了一个自动化流程,用于从CopernicusOpenAccessHub下载Sentinel1卫星图像,使用Sentinelsat库和Snappy(遥感软件SNAP的Python绑定)进行预处理。预处理步骤包括添加轨道文件、去除边界噪声、热噪声消除、校准、TOPSAR解突发、多视处理、子集选择和地形校正。代码适用于特定的日期范围和地理区域,使用了GeoJSON足迹定义。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

一、前言

博客的目标是下载 Sentinel 1 图像,预处理并存储它们以供进一步分析。 由于很多这些分析涉及在各种其他软件之间传播的试验和错误,我试图采用最佳实践并避免复杂性。

首先,您应该在 Copernicus Open Access Hub 上免费注册以下载图像。 Hub 允许我们使用他们的 API 下载图片,但目前每个帐户只能下载两张图片。 为此,我还使用了一个名为 Sentinelsat 的 Python 库,因为它更易于使用。 然后我使用了 Snappy,它是一个非常有用且免费的遥感软件 SNAP 的 python 绑定。 SNAP 有 32 位和 64 位两种版本。

将 SNAP 与现有的 Python 安装绑定需要一些功夫。 由于我已经通过 ArcGIS Pro 安装了 64 位 Python,因此我克隆了一个单独的 Python 环境并在其中安装了 Snappy。 Snappy 帮助我们以编程方式访问 SNAP 的大部分功能。

二、代码

import os
import datetime
import gc
import glob
import snappy
from sentinelsat import SentinelAPI, geojson_to_wkt, read_geojson
from snappy import ProductIO

class sentinel1_download_preprocess():
    def __init__(self, input_dir, date_1, date_2, query_style, footprint, lat=24.84, lon=90.43, download=False):
        self.input_dir = input_dir
        self.date_start = datetime.datetime.strptime(date_1, "%d%b%Y")
        self.date_end = datetime.datetime.strptime(date_2, "%d%b%Y")
        self.query_style = query_style
        self.footprint = geojson_to_wkt(read_geojson(footprint))
        self.lat = lat
        self.lon = lon
        self.download = download

        # configurations
        self.api = SentinelAPI('scihub_username', 'scihub_passwd', 'https://scihub.copernicus.eu/dhus')
        self.producttype = 'GRD'  # SLC, GRD, OCN
        self.orbitdirection = 'ASCENDING'  # ASCENDING, DESCENDING
        self.sensoroperationalmode = 'IW'  # SM, IW, EW, WV

    def sentinel1_download(self):
        global download_candidate
        if self.query_style == 'coordinate':
            download_candidate = self.api.query('POINT({0} {1})'.format(self.lon, self.lat),
                                                date=(self.date_start, self.date_end),
                                                producttype=self.producttype,
                                                orbitdirection=self.orbitdirection,
                                                sensoroperationalmode=self.sensoroperationalmode)
        elif self.query_style == 'footprint':
            download_candidate = self.api.query(self.footprint,
                                                date=(self.date_start, self.date_end),
                                                producttype=self.producttype,
                                                orbitdirection=self.orbitdirection,
                                                sensoroperationalmode=self.sensoroperationalmode)
        else:
            print("Define query attribute")

        title_found_sum = 0
        for key, value in download_candidate.items():
            for k, v in value.items():
                if k == 'title':
                    title_info = v
                    title_found_sum += 1
                elif k == 'size':
                    print("title: " + title_info + " | " + v)
        print("Total found " + str(title_found_sum) +
              " title of " + str(self.api.get_products_size(download_candidate)) + " GB")

        os.chdir(self.input_dir)
        if self.download:
            if glob.glob(input_dir + "*.zip") not in [value for value in download_candidate.items()]:
                self.api.download_all(download_candidate)
                print("Nothing to download")
        else:
            print("Escaping download")
        # proceed processing after download is complete
        self.sentinel1_preprocess()

    def sentinel1_preprocess(self):
        # Get snappy Operators
        snappy.GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
        # HashMap Key-Value pairs
        HashMap = snappy.jpy.get_type('java.util.HashMap')

        for folder in glob.glob(self.input_dir + "\*"):
            gc.enable()
            if folder.endswith(".zip"):
                timestamp = folder.split("_")[5]
                sentinel_image = ProductIO.readProduct(folder)
                if self.date_start <= datetime.datetime.strptime(timestamp[:8], "%Y%m%d") <= self.date_end:
                    # add orbit file
                    self.sentinel1_preprocess_orbit_file(timestamp, sentinel_image, HashMap)
                    # remove border noise
                    self.sentinel1_preprocess_border_noise(timestamp, HashMap)
                    # remove thermal noise
                    self.sentinel1_preprocess_thermal_noise_removal(timestamp, HashMap)
                    # calibrate image to output to Sigma and dB
                    self.sentinel1_preprocess_calibration(timestamp, HashMap)
                    # TOPSAR Deburst for SLC images
                    if self.producttype == 'SLC':
                        self.sentinel1_preprocess_topsar_deburst_SLC(timestamp, HashMap)
                    # multilook
                    self.sentinel1_preprocess_multilook(timestamp, HashMap)
                    # subset using a WKT of the study area
                    self.sentinel1_preprocess_subset(timestamp, HashMap)
                    # finally terrain correction, can use local data but went for the default 
                    self.sentinel1_preprocess_terrain_correction(timestamp, HashMap)
                    # break # try this if you want to check the result one by one
            
    def sentinel1_preprocess_orbit_file(self, timestamp, sentinel_image, HashMap):
        start_time_processing = datetime.datetime.now()
        orb = self.input_dir + "\\orb_" + timestamp

        if not os.path.isfile(orb + ".dim"):
            parameters = HashMap()
            orbit_param = snappy.GPF.createProduct("Apply-Orbit-File", parameters, sentinel_image)
            ProductIO.writeProduct(orbit_param, orb, 'BEAM-DIMAP')  # BEAM-DIMAP, GeoTIFF-BigTiff
            print("orbit file added: " + orb +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + orb)

    def sentinel1_preprocess_border_noise(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        border = self.input_dir + "\\bordr_" + timestamp

        if not os.path.isfile(border + ".dim"):
            parameters = HashMap()
            border_param = snappy.GPF.createProduct("Remove-GRD-Border-Noise", parameters,
                                                    ProductIO.readProduct(self.input_dir +
                                                                          "\\orb_" + timestamp + ".dim"))
            ProductIO.writeProduct(border_param, border, 'BEAM-DIMAP')
            print("border noise removed: " + border +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + border)

    def sentinel1_preprocess_thermal_noise_removal(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        thrm = self.input_dir + "\\thrm_" + timestamp

        if not os.path.isfile(thrm + ".dim"):
            parameters = HashMap()
            thrm_param = snappy.GPF.createProduct("ThermalNoiseRemoval", parameters,
                                                  ProductIO.readProduct(self.input_dir + "\\bordr_" +
                                                                        timestamp + ".dim"))
            ProductIO.writeProduct(thrm_param, thrm, 'BEAM-DIMAP')
            print("thermal noise removed: " + thrm +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + thrm)

    def sentinel1_preprocess_calibration(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        calib = self.input_dir + "\\calib_" + timestamp

        if not os.path.isfile(calib + ".dim"):
            parameters = HashMap()
            parameters.put('outputSigmaBand', True)
            parameters.put('outputImageScaleInDb', False)
            calib_param = snappy.GPF.createProduct("Calibration", parameters,
                                                   ProductIO.readProduct(self.input_dir + "\\thrm_" +
                                                                         timestamp + ".dim"))
            ProductIO.writeProduct(calib_param, calib, 'BEAM-DIMAP')
            print("calibration complete: " + calib +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + calib)

    def sentinel1_preprocess_topsar_deburst_SLC(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        deburst = self.input_dir + "\\dburs_" + timestamp

        if not os.path.isfile(deburst):
            parameters = HashMap()
            parameters.put('outputSigmaBand', True)
            parameters.put('outputImageScaleInDb', False)
            deburst_param = snappy.GPF.createProduct("TOPSAR-Deburst", parameters,
                                                     ProductIO.readProduct(self.input_dir + "\\calib_" +
                                                                           timestamp + ".dim"))
            ProductIO.writeProduct(deburst_param, deburst, 'BEAM-DIMAP')
            print("deburst complete: " + deburst +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + deburst)

    def sentinel1_preprocess_multilook(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        multi = self.input_dir + "\\multi_" + timestamp

        if not os.path.isfile(multi + ".dim"):
            parameters = HashMap()
            parameters.put('outputSigmaBand', True)
            parameters.put('outputImageScaleInDb', False)
            multi_param = snappy.GPF.createProduct("Multilook", parameters,
                                                   ProductIO.readProduct(self.input_dir + "\\calib_" +
                                                                         timestamp + ".dim"))
            ProductIO.writeProduct(multi_param, multi, 'BEAM-DIMAP')
            print("multilook complete: " + multi +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + multi)

    def sentinel1_preprocess_subset(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        subset = self.input_dir + "\\subset_" + timestamp

        if not os.path.isfile(subset + ".dim"):
            WKTReader = snappy.jpy.get_type('com.vividsolutions.jts.io.WKTReader')
            
            # converting shapefile to GEOJSON and WKT is easy with any free online tool
            wkt = "POLYGON((92.330290184197 20.5906091141114,89.1246637610338 21.6316051481971," \
                  "89.0330319081811 21.7802436586492,88.0086282580443 24.6678836192818,88.0857830091018 " \
                  "25.9156771178278,88.1771488779853 26.1480664053835,88.3759125970998 26.5942658997298," \
                  "88.3876586919721 26.6120432770312,88.4105534167129 26.6345128356038,89.6787084683935 " \
                  "26.2383305017275,92.348481691233 25.073636976939,92.4252199249342 25.0296592837972," \
                  "92.487261172615 24.9472465376954,92.4967290851295 24.902213855393,92.6799861774377 " \
                  "21.2972058618174,92.6799346581579 21.2853347419811,92.330290184197 20.5906091141114))"

            geom = WKTReader().read(wkt)
            parameters = HashMap()
            parameters.put('geoRegion', geom)
            subset_param = snappy.GPF.createProduct("Subset", parameters,
                                                    ProductIO.readProduct(self.input_dir + "\\multi_" +
                                                                          timestamp + ".dim"))
            ProductIO.writeProduct(subset_param, subset, 'BEAM-DIMAP')
            print("subset complete: " + subset +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + subset)

    def sentinel1_preprocess_terrain_correction(self, timestamp, HashMap):
        start_time_processing = datetime.datetime.now()
        terr = self.input_dir + "\\terr_" + timestamp

        if not os.path.isfile(terr + ".dim"):
            parameters = HashMap()
            # parameters.put('demResamplingMethod', 'NEAREST_NEIGHBOUR')
            # parameters.put('imgResamplingMethod', 'NEAREST_NEIGHBOUR')
            # parameters.put('pixelSpacingInMeter', 10.0)
            terr_param = snappy.GPF.createProduct("Terrain-Correction", parameters,
                                                  ProductIO.readProduct(self.input_dir + "\\subset_" +
                                                                        timestamp + ".dim"))
            ProductIO.writeProduct(terr_param, terr, 'BEAM-DIMAP')
            print("terrain corrected: " + terr +
                  " | took: " + str(datetime.datetime.now() - start_time_processing).split('.', 2)[0])
        else:
            print("file exists - " + terr)

input_dir = "path_to_project_folder\Sentinel_1"
start_date = '01Mar2019'
end_date = '10Mar2019'
query_style = 'footprint' # 'footprint' to use a GEOJSON, 'coordinate' to use a lat-lon 
footprint = 'path_to_project_folder\bd_bbox.geojson'
lat = 26.23
lon = 88.56

sar = sentinel1_download_preprocess(input_dir, start_date, end_date, query_style, footprint, lat, lon, True) 
# proceed to download by setting 'True', default is 'False'
sar.sentinel1_download()

这段代码可以做很多事情! geojson 文件是使用 ArcGIS Pro 从孟加拉国的一个非常通用的 shapefile 创建的。 有很多免费的在线工具可以将 shapefile 转换为 geojson 和 WKT。

这里使用的处理 Sentinel-1 原始文件的步骤不是最通用的方法,请注意,没有真正的方法。 由于不同的研究需要不同的步骤来准备原始数据,因此您需要遵循自己的步骤。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

倾城一少

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值