Commit 9dfb8b01 authored by Lozac'h Loic's avatar Lozac'h Loic
Browse files

ADD:S1Calibration.py

parent 88609d3c
#!/usr/bin/python
# To change this license header, choose License Headers in Project Properties.
# To change this template file, choose Tools | Templates
# and open the template in the editor.
import os
import xml.etree.ElementTree as ET
import argparse
import otbApplication, math
from osgeo import gdal, osr, ogr
from subprocess import Popen, PIPE
globconfig = {
"urlsrmt":"http://step.esa.int/auxdata/dem/SRTMGL1/",
"srtmhgtzip":"/media/sf_D_DRIVE/partage_vbox/Donnees/SRTMHGTZIP",
"demtif":"/media/sf_D_DRIVE/partage_vbox/Donnees/SRTMHGTZIP/TIF",
"geoid":"/media/sf_D_DRIVE/partage_vbox/Donnees/egm96.grd"
}
namespaces = {
'gml': "http://www.opengis.net/gml"
}
def search_files(directory='.', resolution='S1', extension='SAFE', fictype='d'):
images=[]
extension = extension.lower()
resolution = resolution.lower()
for dirpath, dirnames, files in os.walk(directory):
if fictype == 'f':
for name in files:
# print(os.path.join(dirpath, name) + " test")
if extension and name.lower().endswith(extension) and name.lower().find(resolution) >= 0 :
# print(os.path.join(dirpath, name) + " OK")
abspath = os.path.abspath(os.path.join(dirpath, name))
images.append(abspath)
elif fictype == 'd':
for dirname in dirnames:
# print(os.path.join(dirpath, name) + " test")
if extension and dirname.lower().endswith(extension) and dirname.lower().find(resolution) >= 0 :
# print(os.path.join(dirpath, name) + " OK")
abspath = os.path.abspath(os.path.join(dirpath, dirname))
images.append(abspath)
else:
print("search_files type error")
exit()
return images
def process_command(cmd):
print("Starting : "+" ".join(cmd))
p = Popen(cmd, stdout=PIPE)
# p.wait()
output = p.communicate()[0]
if p.returncode != 0:
print("process failed %d : %s" % (p.returncode, output))
print("#################################################")
return p.returncode
def prepare_srtm_hgt(sarsafe):
#get tiles
tree = ET.parse(os.path.join(sarsafe,"manifest.safe"))
root = tree.getroot()
footprintnode = root.find(".//gml:coordinates",namespaces)
footprintcoordsp = footprintnode.text.split(" ")
footprintlats = []
footprintlons = []
for coord in footprintcoordsp:
coordsp = coord.split(",")
footprintlats.append(float(coordsp[0]))
footprintlons.append(float(coordsp[1]))
llx, lly = min(footprintlons),min(footprintlats)
urx, ury = max(footprintlons),max(footprintlats)
pointLL, pointUR = ogr.Geometry(ogr.wkbPoint), ogr.Geometry(ogr.wkbPoint)
pointLL.AddPoint(llx, lly)
pointUR.AddPoint(urx, ury)
tilesname=[]
print("LowerLeft="+pointLL.ExportToWkt()+";UpperRight="+pointUR.ExportToWkt())
for o in range(math.ceil(pointUR.GetX()) - math.floor(pointLL.GetX())):
for a in range(math.ceil(pointUR.GetY()) - math.floor(pointLL.GetY())):
lat = math.floor(pointLL.GetY()) + a
lon = math.floor(pointLL.GetX()) + o
if lat >= 0 :
hem = 'N'
else :
hem = 'S'
lat = abs(lat)
if lon >= 0 :
grw = 'E'
else :
grw = 'W'
lon = abs(lon)
tilesname.append(hem+f'{lat:02}'+grw+f'{lon:003}')
print("Checking SRTM tiles: "+str(tilesname))
tilefiles = []
if len(tilesname) == 0:
print("bug")
exit()
srtmsuf = ".SRTMGL1.hgt.zip"
for tile in tilesname :
tilef = os.path.join(globconfig["srtmhgtzip"],tile + srtmsuf)
if not os.path.exists(tilef):
print("Download SRTM tile : "+tilef)
cmd=["wget","-P", globconfig['srtmhgtzip'],globconfig["urlsrmt"]+tile + srtmsuf]
rcode = process_command(cmd)
if rcode != 0:
continue
hgtdir = os.path.join(globconfig["srtmhgtzip"],"HGT")
tifdir = os.path.join(globconfig["srtmhgtzip"],"TIF")
tilehgt = os.path.join(hgtdir,tile+".hgt")
tiletif = os.path.join(tifdir,tile+".tif")
if not os.path.exists(tilehgt):
print("Unzip SRTM tile : "+tilehgt)
cmd=['unzip','-d',hgtdir,tilef]
process_command(cmd)
if not os.path.exists(tiletif):
print("Convert SRTM tile : "+tiletif)
cmd=['gdal_translate',tilehgt,tiletif]
process_command(cmd)
print("Ok for SRTM tiles: "+str(tilesname))
def create_incidenceAngle_raster(annofiles,img):
outdir = os.path.dirname(img)
tempdir = os.path.join(os.path.dirname(img),"auxfiles")
imgbs = os.path.basename(img)
imgsp = os.path.basename(img).split("_")
#todo: mettre un dosiier temp
outinc = os.path.join(outdir,"_".join(imgsp[:-1])+"_THETA.TIF")
vrtfile = os.path.join(tempdir,"_".join(imgsp[:-1])+"_THETA.vrt")
csvfile = os.path.join(tempdir,"_".join(imgsp[:-1])+"_THETA.csv")
layername = "_".join(imgsp[:-1])+"_THETA"
if not os.path.exists(tempdir):
os.mkdir(tempdir)
raster = gdal.Open(img)
proj = osr.SpatialReference(wkt=raster.GetProjection())
upx, xres, xskew, upy, yskew, yres = raster.GetGeoTransform()
cols = raster.RasterXSize
rows = raster.RasterYSize
ulx = upx + 0*xres + 0*xskew
uly = upy + 0*yskew + 0*yres
llx = upx + 0*xres + rows*xskew
lly = upy + 0*yskew + rows*yres
lrx = upx + cols*xres + rows*xskew
lry = upy + cols*yskew + rows*yres
urx = upx + cols*xres + 0*xskew
ury = upy + cols*yskew + 0*yres
pointLL, pointUR = ogr.Geometry(ogr.wkbPoint), ogr.Geometry(ogr.wkbPoint)
pointLL.AddPoint(llx, lly)
pointUR.AddPoint(urx, ury)
if not proj.IsGeographic() :
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)
coordTransform = osr.CoordinateTransformation(proj, outSpatialRef)
pointLL.Transform(coordTransform)
pointUR.Transform(coordTransform)
for annofile in annofiles:
tree = ET.parse(annofile)
root = tree.getroot()
geoloclist = root.find(".//geolocationGridPointList")
if geoloclist == None :
print("Can't find geolocationGridPointList tag in xml file.")
exit()
header = True
if os.path.exists(csvfile):
header = False
if not os.path.exists(vrtfile):
with open(vrtfile,'w') as vrt:
vrt.write("<OGRVRTDataSource>\n")
vrt.write(' <OGRVRTLayer name="'+layername+'">\n')
vrt.write(' <SrcDataSource>'+csvfile+'</SrcDataSource>\n')
vrt.write(" <GeometryType>wkbPoint</GeometryType>\n")
vrt.write(" <LayerSRS>WGS84</LayerSRS>\n")
vrt.write(' <GeometryField encoding="PointFromColumns" x="Easting" y="Northing" z="incidence"/>\n')
vrt.write(" </OGRVRTLayer>\n")
vrt.write("</OGRVRTDataSource>\n")
with open(csvfile,'a') as csvf:
if header:
csvf.write("Easting,Northing,incidence\n")
for child in geoloclist:
lat = float(child.find(".//latitude").text)
lon = float(child.find(".//longitude").text)
inc = float(child.find(".//incidenceAngle").text)
csvf.write(str(lon)+","+str(lat)+","+str(inc)+"\n")
cmd = ["gdal_grid","-a","linear","-txe", str(pointLL.GetX()), str(pointUR.GetX()), "-tye", str(pointLL.GetY()), str(pointUR.GetY()), "-outsize", str(cols), str(rows), "-of", "GTiff", "-ot", "Float32", "-l", layername, vrtfile, outinc]
process_command(cmd)
def cal_ortho_sub_mosa_pipeline(sarfiles, ref, outfile):
print("processing.......................................................")
# print("image1 = "+mos1)
# print("image2 = "+mos2)
# print("output = "+out)
apps=[]
i=0
for sarfile in sarfiles:
apps.append(otbApplication.Registry.CreateApplication("SARCalibration"))
apps[i].SetParameterString("in",sarfile)
apps[i].SetParameterString("out", str(i)+"temp.tif")
apps[i].Execute()
apps.append(otbApplication.Registry.CreateApplication("BandMath"))
apps[i+1].AddImageToParameterInputImageList("il",apps[i].GetParameterOutputImage("out"))
apps[i+1].SetParameterString("out", str(i+1)+"temp.tif")
apps[i+1].SetParameterString("exp", "im1b1<=0?0:im1b1")
apps[i+1].Execute()
apps.append(otbApplication.Registry.CreateApplication("OrthoRectification"))
apps[i+2].SetParameterInputImage("io.in",apps[i+1].GetParameterOutputImage("out"))
apps[i+2].SetParameterString("io.out",str(i+2)+"temp.tif")
apps[i+2].SetParameterFloat("opt.gridspacing",40)
apps[i+2].SetParameterString("elev.dem", globconfig["demtif"])
apps[i+2].SetParameterString("elev.geoid", globconfig["geoid"])
apps[i+2].Execute()
apps.append(otbApplication.Registry.CreateApplication("Superimpose"))
apps[i+3].SetParameterInputImage("inm",apps[i+2].GetParameterOutputImage("io.out"))
apps[i+3].SetParameterString("out", str(i+3)+"temp.tif")
apps[i+3].SetParameterString("interpolator","nn")
apps[i+3].SetParameterString("inr",ref)
apps[i+3].Execute()
i+=4
if i == 8 :
apps.append(otbApplication.Registry.CreateApplication("Mosaic"))
apps[i].AddImageToParameterInputImageList("il",apps[3].GetParameterOutputImage("out"))
apps[i].AddImageToParameterInputImageList("il", apps[7].GetParameterOutputImage("out"))
apps[i].SetParameterString("out", str(i)+"temp.tif")
apps[i].Execute()
apps.append(otbApplication.Registry.CreateApplication("ManageNoData"))
apps[i+1].SetParameterInputImage("in", apps[8].GetParameterOutputImage("out"))
apps[i+1].SetParameterString("out", outfile+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
apps[i+1].SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
apps[i+1].SetParameterString("mode", "changevalue")
apps[i+1].ExecuteAndWriteOutput()
else:
apps.append(otbApplication.Registry.CreateApplication("ManageNoData"))
apps[i].SetParameterInputImage("in", apps[3].GetParameterOutputImage("out"))
apps[i].SetParameterString("out", outfile+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
apps[i].SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
apps[i].SetParameterString("mode", "changevalue")
apps[i].ExecuteAndWriteOutput()
print("done.............................................................")
if __name__ == "__main__":
# Make parser object
parser = argparse.ArgumentParser(description=
"""
Compute cloud-masked NDVI from Sentinel2 or Landsat8.
""")
parser.add_argument('-indir', action='store', required=True, help='Directory containing Sentinel1 SNAP Calibrated or Landsat8 NDVI')
parser.add_argument('-inref', action='store', required=True, help='Image references Sentinel-2 Tile, example any S2 NDVI')
parser.add_argument('-zone', action='store', required=True, help='Geographic zone reference for output files naming')
parser.add_argument('-outdir', action='store', required=True, help='Output directory')
# parser.add_argument('-user', action='store', required=True, help='SciHub Copernicus user name')
# parser.add_argument('-psswd', action='store', required=True, help='SciHub Copernicus user password')
args=parser.parse_args()
if not os.path.exists(args.outdir):
os.mkdir(args.outdir)
indir=[]
indir=search_files(args.indir, 'S1', 'SAFE', 'd')
dejafait=[]
for file in indir:
if file in dejafait:
continue
prepare_srtm_hgt(file)
bname = os.path.basename(file).split(".")[0]
splitbname = bname.split("_")
plateforme = splitbname[0]
acqdate = splitbname[4].split("T")[0]
acqtime = splitbname[4].split("T")[1]
outfilebase = os.path.join(args.outdir,"_".join(splitbname[:3])+"_"+args.zone.upper()+"_"+acqdate+"T"+acqtime)
l = [indir.index(i) for i in indir if (acqdate in i) and (plateforme in i)]
if len(l) == 2 :
in1vv = search_files(indir[l[0]], '-vv-', 'tiff', 'f')
in2vv = search_files(indir[l[1]], '-vv-', 'tiff', 'f')
if len(in1vv) !=1 or len(in2vv) !=1:
print("Error: Can't find tiff files!")
continue
sarfiles = [in1vv[0],in2vv[0]]
outfile = outfilebase +"_VV.TIF"
dejafait.append(indir[l[0]])
dejafait.append(indir[l[1]])
if not os.path.exists(outfile) :
cal_ortho_sub_mosa_pipeline(sarfiles, args.inref, outfile)
anno1vv = search_files(indir[l[0]], '-vv-', 'xml', 'f')
anno2vv = search_files(indir[l[1]], '-vv-', 'xml', 'f')
annofiles=[]
for a in anno1vv:
if a.find("calibration") >= 0 or a.find("noise") >= 0:
continue
annofiles.append(a)
for a in anno2vv:
if a.find("calibration") >= 0 or a.find("noise") >= 0:
continue
annofiles.append(a)
if len(annofiles) != 2:
print("Error : Can't find annotation files")
exit()
elif len(l) == 1 :
in1vv = search_files(file, '-vv-', 'tiff', 'f')
if len(in1vv) !=1:
print("Error: Can't find tiff files!")
continue
outfile = outfilebase +"_VV.TIF"
if not os.path.exists(outfile) :
cal_ortho_sub_mosa_pipeline(in1vv, args.inref, outfile)
anno1vv = search_files(file, '-vv-', 'xml', 'f')
annofiles=[]
for a in anno1vv:
if a.find("calibration") >= 0 or a.find("noise") >= 0:
continue
annofiles.append(a)
if len(annofiles) != 1:
print("Error : Can't find annotation files")
exit()
else :
print("Error on index :"+l)
create_incidenceAngle_raster(annofiles,outfile)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment