更新2023.12.15
更新
如果要使用这个进行批处理可以注意以下事情:
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