Commit eb714a5c authored by Lozac'h Loic's avatar Lozac'h Loic

Update: ajout de la creation d'un mnt et du slope

parent 199ee425
......@@ -58,23 +58,44 @@ def process_command(cmd):
print("process failed %d : %s" % (p.returncode, output))
print("#################################################")
return p.returncode
def get_img_extend(img):
raster = gdal.Open(img)
proj = osr.SpatialReference(wkt=raster.GetProjection())
def prepare_srtm_hgt(sarsafe):
#get tiles
tree = ET.parse(os.path.join(sarsafe,"manifest.safe"))
root = tree.getroot()
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
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]))
pointLL, pointUR = ogr.Geometry(ogr.wkbPoint), ogr.Geometry(ogr.wkbPoint)
pointLL.AddPoint(llx, lly)
pointUR.AddPoint(urx, ury)
if not proj.IsGeographic() :
llx, lly = min(footprintlons),min(footprintlats)
urx, ury = max(footprintlons),max(footprintlats)
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)
coordTransform = osr.CoordinateTransformation(proj, outSpatialRef)
pointLL.Transform(coordTransform)
pointUR.Transform(coordTransform)
return pointLL, pointUR, cols, rows
def find_srtm_hgt_name(llx, lly , urx, ury ):
pointLL, pointUR = ogr.Geometry(ogr.wkbPoint), ogr.Geometry(ogr.wkbPoint)
pointLL.AddPoint(llx, lly)
......@@ -98,6 +119,27 @@ def prepare_srtm_hgt(sarsafe):
grw = 'W'
lon = abs(lon)
tilesname.append(hem+f'{lat:02}'+grw+f'{lon:003}')
return tilesname
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)
tilesname = find_srtm_hgt_name(llx, lly , urx, ury )
print("Checking SRTM tiles: "+str(tilesname))
tilefiles = []
......@@ -107,7 +149,7 @@ def prepare_srtm_hgt(sarsafe):
srtmsuf = ".SRTMGL1.hgt.zip"
tiftiles=[]
for tile in tilesname :
tilef = os.path.join(globconfig["srtmhgtzip"],tile + srtmsuf)
if not os.path.exists(tilef):
......@@ -133,18 +175,54 @@ def prepare_srtm_hgt(sarsafe):
cmd=['gdal_translate',tilehgt,tiletif]
process_command(cmd)
tiftiles.append(tiletif)
print("Ok for SRTM tiles: "+str(tilesname))
return tiftiles
def create_zone_dem_and_slope(ref,zone,outdir):
#get tiles
pointLL, pointUR, cols, rows = get_img_extend(ref)
tilesname = find_srtm_hgt_name(pointLL.GetX(), pointLL.GetY() , pointUR.GetX(), pointUR.GetY() )
tilefiles = []
for tile in tilesname:
tifdir = os.path.join(globconfig["srtmhgtzip"],"TIF")
tilefiles.append(os.path.join(tifdir,tile+".tif"))
print("Generating DEM file on ref:"+ref)
outdem = os.path.join(outdir,"SRTM_"+zone+".TIF")
mosapp = otbApplication.Registry.CreateApplication("Mosaic")
mosparams = {"il":tilefiles, "out":"mosatemp.tif"}
mosapp.SetParameters(mosparams)
mosapp.Execute()
supapp = otbApplication.Registry.CreateApplication("Superimpose")
supapp.SetParameterString("inr", ref)
supapp.SetParameterInputImage("inm", mosapp.GetParameterOutputImage("out"))
supapp.SetParameterString("out", outdem)
supapp.ExecuteAndWriteOutput()
#Compute slope
print("Generating slope file...")
outslope = os.path.join(outdir,"SLOPE_"+zone+".TIF")
cmd=['gdaldem', 'slope', '-p', outdem, outslope]
process_command(cmd)
return outdem
def create_incidenceAngle_raster(annofiles,img):
def create_incidenceAngle_raster(annofiles,dem, 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")
outinc = os.path.join(tempdir,"_".join(imgsp[:-1])+"_THETASURF.TIF")
outincfromdem = 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"
......@@ -152,37 +230,7 @@ def create_incidenceAngle_raster(annofiles,img):
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)
pointLL, pointUR, cols, rows = get_img_extend(img)
for annofile in annofiles:
......@@ -214,19 +262,40 @@ def create_incidenceAngle_raster(annofiles,img):
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")
lat = child.find(".//latitude").text
lon = child.find(".//longitude").text
alti = float(child.find(".//height").text)
incdem = float(child.find(".//incidenceAngle").text)
inc = incdem - alti/693000.
csvf.write(lon+","+lat+","+str(incdem)+"\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)
app41 = otbApplication.Registry.CreateApplication("Superimpose")
# The following lines set all the application parameters:
app41.SetParameterString("inm", outinc)
app41.SetParameterString("inr", dem)
app41.SetParameterString("interpolator","nn")
app41.SetParameterString("out", "temp41.tif")
print("Launching... Resampling")
# The following line execute the application
app41.Execute()
print("End of Resampling \n")
appS = otbApplication.Registry.CreateApplication("BandMath")
appS.SetParameterStringList("il", [dem])
appS.AddImageToParameterInputImageList("il", app41.GetParameterOutputImage("out"))
# Define Input im2: Band Red (B4)
appS.SetParameterString("out", outincfromdem)
appS.SetParameterString("exp", "im2b1 - im1b1/693000" )
appS.ExecuteAndWriteOutput()
def cal_ortho_sub_mosa_pipeline(sarfiles, ref, outfile):
......@@ -309,10 +378,16 @@ if __name__ == "__main__":
if not os.path.exists(args.outdir):
os.mkdir(args.outdir)
demdir = os.path.join(args.outdir,"DEM")
if not os.path.exists(demdir):
os.mkdir(demdir)
indir=[]
indir=search_files(args.indir, 'S1', 'SAFE', 'd')
zonedemfile = create_zone_dem_and_slope(args.inref,args.zone,demdir)
dejafait=[]
for file in indir:
if file in dejafait:
......@@ -395,4 +470,7 @@ if __name__ == "__main__":
print("Error on index :"+l)
create_incidenceAngle_raster(annofiles,outfile)
create_incidenceAngle_raster(annofiles,zonedemfile,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