Commit 576966d7 authored by Lozac'h Loic's avatar Lozac'h Loic
Browse files

Delete batch_create_mv_metadata.py

parent 5cdd7237
#!/usr/bin/python
__author__ = "Loic Lozach"
__date__ = "$Dec 14, 2018 12:40:21 PM$"
import json, os, sys, argparse
from subprocess import Popen, PIPE
try:
from osgeo import ogr, osr, gdal
except:
sys.exit('ERROR: cannot find GDAL/OGR modules')
import otbApplication
def create_dict_assets(smhref, qklhref):
mvassets = {
"soil_moisture":{
"href":smhref,
"type":"image/vnd.stac.geotiff"
}, "thumbnail":{
"href":qklhref,
"type":"image/png"
}
}
return mvassets
def create_dict_properties(datetime,platform,gsd,epsg,title,description):
#gsd => resolution en metre
mvproperties = {
"datetime":datetime,
"title":title,
"description":description,
"platform":platform,
"eo:gsd":gsd,
"eo:epsg":epsg
}
return mvproperties
def create_dict_geometry(coordinates,geomtype="Polygon"):
mvgeom = {
"type":geomtype,
"coordinates":[coordinates]
}
return mvgeom
def create_dict_image(imgid, geometry, properties, assets, imgtype="Feature"):
mvdata = {
"id":imgid,
"type":imgtype,
"geometry":geometry,
"properties":properties,
"assets":assets
}
return mvdata
def search_files(directory='.', extension='tif', resolution='MV_'):
images=[]
extension = extension.lower()
resolution = resolution.lower()
for dirpath, dirnames, files in os.walk(directory):
for name in files:
if dirpath.upper().find("_PROD") >= 0 or dirpath.upper().find("_TMP") >= 0 :
continue
if extension and name.lower().endswith(extension) and name.lower().find(resolution) >= 0 :
upname = name.upper()
if upname != name :
upabspath = os.path.abspath(os.path.join(dirpath, upname))
abspath = os.path.abspath(os.path.join(dirpath, name))
os.rename(abspath, upabspath)
name = upname
#print(os.path.join(dirpath, name))
abspath = os.path.abspath(os.path.join(dirpath, name))
images.append(abspath)
return images
def generate_quicklook(img, lut, quick):
# The following line creates an instance of the Superimpose application
app1 = otbApplication.Registry.CreateApplication("ColorMapping")
# The following lines set all the application parameters:
app1.SetParameterString("in", img)
app1.SetParameterString("out", "temp1.tif")
app1.SetParameterString("op", "labeltocolor")
app1.SetParameterString("method", "custom")
app1.SetParameterString("method.custom.lut", lut)
print("Launching... Resampling")
# The following line execute the application
app1.Execute() #ExecuteAndWriteOutput()
print("End of Resampling \n")
app2 = otbApplication.Registry.CreateApplication("Quicklook")
# The following lines set all the application parameters:
app2.SetParameterInputImage("in", app1.GetParameterOutputImage("out"))
app2.SetParameterString("out", quick)
app2.SetParameterInt("sx", 1024)
print("Launching... Resampling")
# The following line execute the application
app2.ExecuteAndWriteOutput() #ExecuteAndWriteOutput()
print("End of Resampling \n")
def reprojectRaster(absfile):
temp_file = absfile + ".old"
os.rename(absfile, temp_file)
p = Popen(['gdalwarp', '-t_srs', 'EPSG:4326', temp_file, absfile], stdout=PIPE)
# p.wait()
output = p.communicate()[0]
if p.returncode != 0:
print("gdalwarp failed %d : %s" % (p.returncode, output))
with open("gdalwrap_log.err",'a') as err:
err.write("######################################################################################################")
err.write(absfile)
err.write(output)
err.write("######################################################################################################")
os.remove(absfile)
os.rename(temp_file, absfile)
return 1
print("gdalwarp succeeded on : "+absfile)
os.remove(temp_file)
return 0
if __name__ == "__main__":
# Make parser object
parser = argparse.ArgumentParser(description=
"""
Create geojson metadata recursively for Soil Moisture images collections
""")
#####################################################
driver = gdal.GetDriverByName('GTiff')
driver.Register()
#####################################################
parser.add_argument('-mvdir', action='store', required=True, help='Directory containing Soil Moisture data')
parser.add_argument('--quicklook', dest='quicklook', action='store_true', help='Generate quicklook images (default True)')
parser.add_argument('--no-quicklook', dest='quicklook', action='store_false')
parser.set_defaults(quicklook=True)
parser.add_argument('-lut', action='store', required=False, help='LUT file for quicklook generation')
parser.add_argument('--overwrite', dest='overwrite', action='store_true', help='Overwrite already existing files')
parser.add_argument('--no-overwrite', dest='overwrite', action='store_false', help='Overwrite already existing files (default False)')
parser.set_defaults(overwrite=False)
parser.add_argument('--reproj', dest='reproj', action='store_true', help='Reprojection in WGS84 if not (default True)')
parser.add_argument('--no-reproj', dest='reproj', action='store_false', help='Reprojection in WGS84 if not')
parser.set_defaults(reproj=True)
args=parser.parse_args()
if args.quicklook:
if None == args.lut:
print("LUT is required for quicklook generation")
exit
mvtifs=[]
mvtifs=search_files(args.mvdir)
if len(mvtifs) == 0 :
print("No soil moisture image found")
exit
for file in mvtifs:
print("################################################################")
print("Processing file " + file)
absfile = os.path.abspath(file)
spltpath = absfile.split("/")
geojsonfile = absfile[:-3]+"JSON"
quicklookfile = absfile[:-3]+"PNG"
if os.path.exists(geojsonfile) and not args.overwrite :
print(geojsonfile + " already exists! Skipping file.")
continue
filename = os.path.basename(file)
spltfilename = filename.split("_")
date = spltfilename[len(spltfilename)-1][:-4]
date_in_format = date[0:4]+"-"+date[4:6]+"-"+date[6:8]+date[8:11]+":"+date[11:13]+":"+date[13:15]+"Z"
mv_plateform = spltfilename[1]
ind = spltpath.index("SoilMoisture")
mv_url = "https://products.thisme.cines.teledetection.fr/" + '/'.join(spltpath[ind+1:])
qkl_url = mv_url[:-3]+"PNG"
mv_title=spltpath[ind+1]
mv_description="Sentinel-1 derived soil Moisture product for #"+mv_title+" at Plot scale (S2MP)"
raster = gdal.Open(absfile)
proj = osr.SpatialReference(wkt=raster.GetProjection())
mv_epsg = int(proj.GetAttrValue('AUTHORITY',1))
if args.reproj and mv_epsg != 4326 :
print("Reprojection needed. Launching gdalwrap.")
raster = None
projerr = reprojectRaster(absfile)
if projerr:
continue
raster = gdal.Open(absfile)
proj = osr.SpatialReference(wkt=raster.GetProjection())
mv_epsg = int(proj.GetAttrValue('AUTHORITY',1))
transform = raster.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]
rows = raster.RasterYSize
cols = raster.RasterXSize
lastX = xOrigin + cols *pixelWidth
lastY = yOrigin + rows *pixelHeight
xCenter = xOrigin + cols/2 *pixelWidth
yCenter = yOrigin + rows/2 *pixelHeight
mv_coordinates = [[xOrigin,lastY],[xOrigin,yOrigin],[lastX,yOrigin],[lastX,lastY],[xOrigin,lastY]]
mv_properties = create_dict_properties(date_in_format,mv_plateform,pixelWidth,mv_epsg,mv_title,mv_description)
mv_assets = create_dict_assets(mv_url, qkl_url)
mv_geometry = create_dict_geometry(mv_coordinates)
mv_geojson = create_dict_image(filename[:-3], mv_geometry, mv_properties, mv_assets)
print("Writing json file.")
with open(geojsonfile,'w') as jsonf:
json.dump(mv_geojson, jsonf, indent=4)
if not args.quicklook :
continue
print("Generating quicklook")
generate_quicklook(absfile, args.lut, quicklookfile)
print("Process done.")
print("################################################################")
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