I want to store a large LIDAR dataset as a VRT so I can quickly query for elevation values in the entire state of Massachussets. However, I am having trouble actually creating the mosaic. The naive way of performing this operation, with gdal_merge.py is impractical since the dataset is 21GB of memory. Because of the size of the dataset, it also spans multiple UTM zones so I need to translate at least some of the geotiffs. I wrote a brief python script to facilitate this:
from glob import glob from os import path, system from tqdm import tqdm from osgeo import gdal, osr from gdalconst import * import subprocess #proj4 = '+proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs' #tgt_proj4 = osr.SpatialReference() #tgt_proj4.ImportFromProj4(proj4) ds = gdal.Open("LIDAR/19_03344719.tif", GA_ReadOnly) prj = ds.GetProjection() tgt_srs = osr.SpatialReference(wkt=prj) files = list(glob("LIDAR/*")) files.sort(reverse=True) print(tgt_srs.ExportToWkt()) for filename in tqdm(files):e newfilename = path.join("LIDAR_TRANS", path.splitext(path.basename(filename)) + ".tif") subprocess.call(['gdalwarp', filename, newfilename, '-t_srs', tgt_srs.ExportToProj4(), '-multi', '-q'])
I then simply merge all the items in LIDAR_TRANS into a giant VRT. However, I had issues doing this as I keeping large gaps in my mosaic VRT as shown below
As such, I was wondering how my pipeline needs to be changed to avoid these gaps so I can properly use the dataset?