5

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))[0] + ".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

Image Error

As such, I was wondering how my pipeline needs to be changed to avoid these gaps so I can properly use the dataset?

  • Gdalbuildvrt gdal.org/gdalbuildvrt.html is often a better tool that gdal_merge.py. Try how it behaves. However, virtual rasters with very many source images tend to be slow. – user30184 Jan 4 '17 at 22:00
  • The issue is that they span different SRS and gdalbuildvrt cannot handle different srs. That's currently what I used but transferring them into the proper srs so I could run gdalbuildvrt created these gaps. – Skylion Jan 4 '17 at 22:12
  • 2
    It looks like the black areas are not gaps but they are nodata areas which appear because of warping. Perhaps with improving the gdalwarp parameters you can make nodata transparent. – user30184 Jan 5 '17 at 6:13
  • Link to the MassGIS LiDAR Terrain Portal. Looks like a terrabyte or more of LAS ( or LAZ!? ) so making a single raster from the raw lidar doesn't look like an easy task. – Logan Byers Jan 6 '17 at 15:30
1

I am quite certain that the "gaps" you are seeing are just the overlay of NoData borders atop good data on rasters below them. The sharp spikes of good data in the corners are the smoking guns.

You could set the nodata value for each raster in the warping stage (-dstnodata) or at the VRT build step (-vrtnodata). The warping stage is probably preferred.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.