from osgeo import gdal import os from datetime import datetime as dt def gen_xml(wmts,output='~/output.xml',type_flux='WMTS'): src_ds = '%s:%s' % (type_flux,wmts) gdal.Translate( output, src_ds, format=type_flux ) return output def gen_mbtiles(xml,output,format='mbtiles'): init = dt.now() creation_options = ["TILE_FORMAT=JPEG", "QUALITY=100","MINZOOM=12","MAXZOOM=20"] gdal.SetCacheMax(1080) gdal.Translate( output, xml, width=6459465-6458808, height=914340-913549, # width=750, # height=750, # widthPct=50, # heightPct=50, # projWin=[836250,6534274,965160,6403454], projWin=[913549,6459465,914340,6458808], projWinSRS='EPSG:2154', # outputBounds=[913549,6459465,914340,6458808], outputSRS='EPSG:2154', format=format.upper(), creationOptions=creation_options, callback=gdal.TermProgress ) print(dt.now()-init) def gen_mbtiles2(xml,output,format='mbtiles'): # projwin for test : 913549 6459465 914340 6458808 # projwin Isère : 836250 6534274 965160 6403454 init = dt.now() cmd = ''' gdal_translate -outsize 75000 75000 -projwin 836250 6534274 965160 6403454 -projwin_srs "EPSG:2154" -of %s -r cubic -co "COMPRESS=JPEG" -co "TILE_FORMAT=JPEG" -co "QUALITY=80" -co "MINZOOM=12" -co "MAXZOOM=22" -co "PHOTOMETRIC=YCBCR" --config GDAL_CACHEMAX 512 "%s" "%s" ''' % (format,xml,output) print(cmd) os.system(cmd) print(dt.now()-init) def gdalwarp(vrt,mbtiles): # exemple : Découpage d'un raster à partir d'un polygon # gdalwarp -overwrite -s_srs "EPSG:2154" -t_srs "EPSG:2154" -of GTiff -cutline "PG:dbname='azalee' host=91.134.194.221 port=5432 sslmode=disable user='cgeier' password='adm1n*bdCen'" -cl "ref_territoire.isere_platiere" -crop_to_cutline -multi -co "COMPRESS=JPEG" -co "JPEG_QUALITY=75" "/media/colas/Disk2/5_BDD/ZH_prob/ZH_probS_Auvergne-Rhone-Alpes/ZH_probS_region84.tif" "/media/colas/Disk2/5_BDD/ZH_prob/ZH_probS_Auvergne-Rhone-Alpes/ZH_probS_dept38.tif" cmd = ''' gdalwarp -of MBTiles -s_srs epsg:2154 -t_srs epsg:3857 %s %s ''' % (vrt,mbtiles) print(cmd) os.system(cmd) def gdaladdo(mbtiles): init = dt.now() cmd = ''' gdaladdo -r cubic --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config USE_RRD YES --config GDAL_NUM_THREADS ALL_CPUS %s 2 4 8 16 ''' % mbtiles print(cmd) os.system(cmd) print(dt.now()-init) if __name__ == "__main__": flux = 'https://data.geopf.fr/wmts?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetCapabilities&dpiMode=7&format=image/jpeg&layers=HR.ORTHOIMAGERY.ORTHOPHOTOS&styles=normal&tileMatrixSet=PM,layer=HR.ORTHOIMAGERY.ORTHOPHOTOS' flux = 'https://wxs.ign.fr/ortho/geoportail/wmts?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetCapabilities&dpiMode=7&format=image/jpeg&layers=HR.ORTHOIMAGERY.ORTHOPHOTOS&styles=normal&tileMatrixSet=PM,layer=HR.ORTHOIMAGERY.ORTHOPHOTOS' flux = 'https://tiles.craig.fr/ign?VERSION=1.3.0&crs=EPSG:2154&featureCount=10&format=image/png&layers=scans_ign&maxHeight=256&maxWidth=256' flux = 'http://tiles.craig.fr/ortho/service?VERSION=1.3.0&crs=EPSG:2154&featureCount=10&format=image/jpeg&layers=ortho_2021&maxHeight=256&maxWidth=256' # flux = 'http://mt.google.com/vt/lyrs=s&x=${x}&y=${y}&z=${z}' PATH = '/media/colas/Disk2' xml = 'googlemaps.xml' xml = 'craig.xml' xml = 'output.xml' xml = 'ortho_craig.xml' vrt = 'ortho_craig.vrt' mbtile = 'satellite.mbtiles' mbtile = 'scan_craig.mbtiles' mbtile = 'ign_ortho2024_38.mbtiles' mbtile = 'ortho_craig.mbtiles' out_xml = gen_xml( wmts=flux, output=os.path.join(PATH,xml), type_flux='wms') # gen_mbtiles2(os.path.join(PATH,xml),os.path.join(PATH,mbtile),format='mbtiles') # gen_mbtiles2(os.path.join(PATH,xml),os.path.join(PATH,vrt),format='vrt') # gdalwarp(os.path.join(PATH,vrt),os.path.join(PATH,mbtile)) gdaladdo(os.path.join(PATH,mbtile)) # ds = gdal.Open(os.path.join(PATH,vrt)) # gdal_translate -outsize 50% 50% -projwin 631397 5672590 639669 5659275 -of MBTILES -co "COMPRESS=YES" -co "TILE_FORMAT=JPEG" -co "QUALITY=80" "/media/colas/Disk2/output.xml" "/media/colas/Disk2/ign_ortho2024_38.mbtiles"