python pyroSAR对Sentinel-1数据进行预处理

更新2023.12.15

另外一个完全基于python的库是：GitHub - aalling93/Sentinel_1_python: Sentinel-1 SAR in Python. Work with Metadata, Download Images, Load image, Calibrate the images, Speckle Filter the images. In short: Work with Sentinel 1 with only python...Sentinel-1 SAR in Python. Work with Metadata, Download Images, Load image, Calibrate the images, Speckle Filter the images. In short: Work with Sentinel 1 with only python... - GitHub - aalling93/Sentinel_1_python: Sentinel-1 SAR in Python. Work with Metadata, Download Images, Load image, Calibrate the images, Speckle Filter the images. In short: Work with Sentinel 1 with only python...https://github.com/aalling93/Sentinel_1_python

---

更新

如果要使用这个进行批处理可以注意以下事情：

1、如果程序跑了很久那一景都没有跑出来也没报错，可能是在下载轨道文件（需梯子）非常缓慢，一定要登梯子。

2、如果报错可能是需要SNAP更新，打开SNAP更新就好（最近都需要梯子）。注意更新完需要重新打开SNAP！

3、pyroSAR本身自己的问题也可能导致报错，重新pip install --update这个库一下，更新完之后还是需要按照我下面写的修改，这样才能有正确的地理信息。

---

查看了目前能够对Sentinel-1进行处理的python库，发现除了SNAP自带的snappy和另一个完全独立于SNAP的python库外，还有一个pyroSAR，基于SNAP本身的XML，但是进行了优化，平均处理每一景2-3min (另一台一般的机器4-6min)。目前网上代码很少，因此做一个分享。

geocode函数能够集合了一些常用的预处理步骤，读取下载的zip文件，直接生成入射角、HH、HV等，每个单独一个tif文件。最后存储的文件命名自动形成，不能更改。处理出来的无效值为nodata，在ENVI中可以显现，但是后面利用gdal读取，nodata的值是0.0，因此需要在后续处理中注意。

from datetime import datetime
from pyroSAR.snap.util import geocode,ID,identify,sub_parametrize

start_time = datetime.now()
file=r'Z:\download\Sentinel\2016\1\S1A_EW_GRDM_1SDH_20160101T013526_20160101T013630_009294_00D6C2_EF1A.zip'
processed_path=r'Z:\download\process\2016\1'
target_resolution=40
terrain_flat_bool=False
remove_therm_noise_bool=True
fileid = identify(file)
# s1file=ID(fileid)
corners=fileid.getCorners()
subsetnode=sub_parametrize(fileid,geometry=corners)
# generate sequence of integer coordinates marking the tie points of the overlapping hgt tiles

geocode(infile=file, outdir=processed_path, spacing=int(target_resolution), polarizations=['HH','HV'],refarea='sigma0',t_srs=int(3413), scaling='dB',
        clean_edges=True,
            terrainFlattening=terrain_flat_bool, removeS1ThermalNoise=remove_therm_noise_bool,export_extra=['localIncidenceAngle'],demName='ACE30', nodataValueAtSea=False)

interval_time = datetime.now()
print(" processed in " +str(interval_time - start_time) + " Hr:min:sec")

我发现如果用原始的代码，最后得到的tif没有正确的projection。

真正输出用的是geocode函数里调用的writer函数，因此修改了pyroSAR/snap/auxil.py文件里面writer函数里的gdal_translate这一行，替换为以下

注意，在writer函数中，新增t_srs这一个变量，在geocode调用writer函数时传递进去就行

from osgeo import osr
#这一行需要增加

if t_srs is None:
    gdal_translate(src=item, dst=name_new, **translateoptions)   #这一行是源代码
else:
     # 下面到else全是自己写的
    src_ds = gdal.Open(item)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(t_srs)
    src_ds.SetProjection(srs.ExportToWkt())
    driver = gdal.GetDriverByName('GTiff')
    dst_ds = driver.CreateCopy(name_new, src_ds, strict=0,
                                        options=["TILED=YES", "COMPRESS=PACKBITS"])
    # Once we're done, close properly the dataset
    dst_ds = None
    src_ds = None
    driver=None