#!/usr/bin/python3 __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 = { "tif":{ "href":smhref, "type":"image/vnd.stac.geotiff" }, "thumbnail":{ "href":qklhref, "type":"image/png", "roles":["thumbnail"] } } 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",epsg=4326): mvgeom = { "type":geomtype, "coordinates":[coordinates], "crs":{ "type":"name", "properties":{ "name":"epsg:"+str(epsg) } } } 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, minsx=1024): tmpf = os.path.join(os.path.dirname(quick),"tmpfile.tif") if os.path.exists(tmpf): os.remove(tmpf) # 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", tmpf) app1.SetParameterString("op", "labeltocolor") app1.SetParameterString("method", "custom") app1.SetParameterString("method.custom.lut", lut) print("Launching... Resampling") # The following line execute the application app1.ExecuteAndWriteOutput() #ExecuteAndWriteOutput() print("End of Resampling \n") app2 = otbApplication.Registry.CreateApplication("Quicklook") # The following lines set all the application parameters: app2.SetParameterString("in", tmpf) app2.SetParameterString("out", quick) app2.SetParameterInt("sx", minsx) print("Launching... Resampling") # The following line execute the application app2.ExecuteAndWriteOutput() #ExecuteAndWriteOutput() print("End of Resampling \n") if os.path.exists(tmpf): os.remove(tmpf) 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 THISME image collections """) ##################################################### driver = gdal.GetDriverByName('GTiff') driver.Register() ##################################################### parser.add_argument('-coll', choices=['mv','bf'], default='mv', required=False, help='Choose collection between SoilMoisture (mv), BuildingFootprint (bf), default mv') parser.add_argument('-colldir', action='store', required=True, help='Directory containing collection data') parser.add_argument('-lut', action='store', required=True, help='LUT file for quicklook generation') parser.add_argument('-overwrite', choices=['no', 'all',"json","png"], default='no', required=False, help='[Optional] Overwrite already existing files, default no') args=parser.parse_args() mvtifs=[] if args.coll == "mv": mvtifs=search_files(args.colldir) elif args.coll == "bf": mvtifs=search_files(args.colldir,"tif","BF_") else: print("Error!") exit() if len(mvtifs) == 0 : print("No collection image found") exit for file in mvtifs: cols=None print("################################################################") print("Processing file " + file) absfile = os.path.abspath(file) spltpath = absfile.split("/") geojsonfile = absfile[:-3]+"JSON" quicklookfile = absfile[:-3]+"PNG" print("Processing file " + geojsonfile) if not os.path.exists(geojsonfile) or args.overwrite == "all" or args.overwrite == "json": filename = os.path.basename(file) spltfilename = filename[:-4].split("_") date = spltfilename[3] if date.find('T') > 0 : date_in_format = date[0:4]+"-"+date[4:6]+"-"+date[6:8]+date[8:11]+":"+date[11:13]+":"+date[13:15]+"Z" else : date_in_format = date[0:4]+"-"+date[4:6]+"-"+date[6:8]+"T00:00:00Z" #ind = spltpath.index("_PROD") mv_url = "https://products.thisme.cines.teledetection.fr/" + '/'.join(spltpath[-2:]) qkl_url = mv_url[:-3]+"PNG" if args.coll == "mv": #mv_title=spltpath[ind+1] mv_title=spltfilename[2] mv_description="Sentinel-1 derived soil Moisture product for "+mv_title+" at Plot scale (S2MP)" mv_plateform = spltfilename[1] elif args.coll == "bf": mv_title=spltfilename[2] mv_description="Buildings footprint from Spot-6/7 images for "+mv_title[:-3] mv_plateform = spltfilename[1] else: print("Error!") exit() 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, geomtype="Polygon",epsg=mv_epsg) 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) else: print("Already exists. Passing. ") print("Processing file " + quicklookfile) if not os.path.exists(quicklookfile) or args.overwrite == "all" or args.overwrite == "png" : minsx = 1024 if not cols: raster = gdal.Open(absfile) cols = raster.RasterXSize if cols < 1024 : minsx = cols print("Generating quicklook") generate_quicklook(absfile, args.lut, quicklookfile, minsx=minsx) print("Process done.") print("################################################################")