How to workaround large mosaic process which is failing?

Question

I need to mosaic about 550Gb of tif imagery together and the software I have tried keeps failing. The area has been split into zones so that the smallest has approx 200 tiles.

I have used the latest versions of ERDAS (Imagine and Mapper), ArcINFO and Global Mapper on a 3.30 gigahertz Intel Xeon E31245, DELL, 16GB RAM, 64-bit Win 7 Professional. Mullti-core (4 total), Hyper-threaded (8 total) machine. My C has 700GB free and D has 1.5TB.

I am looking into using Grass (never have before) but i.image.mosaic only seems to only handle 4 files...some of mine have 600 tiles. Any other options or opensource software to try?

Sorry should add that we can't use a mosaic dataset (or equivalent in other software) as we need to create zones with defined no-data areas as ecw's so that they can be opened in any GIS software and combined with lower resolution/older data when new data doesn't exist seamlessly.

An example of how some mosaiced files look in different sofware. Global Mapper/ERDAS are fine but it's not correct in arcgis.

--- OLDER INFO---

Sorry for the rough drawing. So having the colored areas as 5 zones will minimize the no data areas in the larger AOI.

In arcgis the code is as follows (this is run as a model and not in python as I can't get it to take the tifList input).

arcpy.MosaicToNewRaster_management(tifList+";" +mask,RootOutput,"Tile1.tif","PROJCS['GDA_1994_MGA_Zone_55',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',147.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]","16_BIT_UNSIGNED","0.5","3","MAXIMUM","#")

# Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
# The following inputs are layers or table views: "test2"

arcpy.CopyRaster_management(OutputFile,RootOutput+"Tile1b.tif","#","256","256","NONE","NONE","16_BIT_UNSIGNED")

where tifList should be read in from a csv file but this didn't work in python so I am running the above in a model instead...

I have 1.5TB+ free space on my drive but the process crashes with a 9999 error.

Would even 100 tiles process? -i.e should we look at breaking the zones up further?

Answer 1

I will have to 2nd @blah238's suggestions of using some other method of data access than creating a single mosaiced image. A simple guess would say there is not a desktop out there that could handle the amount of data you would have to process in order to mosaic all of those tiles.
To break it down, there are probably two places where you are running out of space.

1. The first is going to be in your RAM buffer. In order to mosaic data, you are merging everything into a single file, which means ideally, that entire file would be held in memory. You don't have 550GB of RAM, meaning there will be read/write from virtual memory. That is enough to crash the process right there.
2. The other issue is likely to be hard disk space. Many of the raster operations create temp files in your "workspace" directory that are quite large. There may be 2 or 3 or more of these before the final dataset is even written, so it is conceivable to use up all your disk space during processing.

Now, for other solutions. As mentioned in the comments above, there is the option to create a Mosaic Dataset. This dataset will allow you to not only treat all of the individual tiles as a single seamless image, but it also maintains metadata about the individual tiles contained within. It also lets you perform raster operations such as Hillshade.

The other option that I would recommend, based on your comment about wanting to have the zones separated would be to create a Raster Catalog. A Raster Catalog is essentially a group layer. You can add multiple raster datasets to it. They can be managed in a geodatabase, and import the rasters, or simply create an unmanaged dataset where the Raster Catalog maintains paths to the original raster datasets. When you load this layer into ArcMap, you can set display properties to only load in a certain number of raster tiles at once, or set the display scale and resolution.
I am currently using a raster catalog to tile a 100+GB set of aerial photos. The performance is very good. If you are looking for a different type of data storage simply for the means of managing a large number of tiles, then I would really recommend it.

Here is code that you could use to create a Raster Catalog and then import a workspace of tiles into it:

import arcpy
import os

newdir = r"c:\data"
dbname = "Aerialphotos.gdb"
gdbdir = os.path.join(newdir, dbname)
rcat = "AerialCatalog"

arcpy.CreateRasterCatalog_management(gdbdir, rcat,
                                     "NAD 1983 StatePlane California VI FIPS 0406 (US Feet).prj", 
                                     "NAD 1983 StatePlane California VI FIPS 0406 (US Feet).prj",
                                     "MAX_FILE_SIZE_4GB", "1000", "3000", "9000",
                                     "UNMANAGED", "")

#Load all raster datasets in workspace to Raster Catalog
rcatdir = os.path.join(gdbdir, rcat)
rastertiledir = os.path.join(newdir, "tiles")

arcpy.WorkspaceToRasterCatalog_management(rastertiledir, rcatdir,
                                          "INCLUDE_SUBDIRECTORIES",
                                          "PROJECT_ONFLY")

Hope this helps!

------------- Edit

Here is a graphic of the tiles handled by my raster catalog. Note that you can choose to have wireframes shown or the raster data. The raster catalog has an attribute table that you can add fields to, for instance if you wanted to add zone designations as in your graphic. Then, you could choose to only show those rasters in a particular zone.
When printing in a graphic from layout view, the full resolution of the rasters is used, so there is no loss of quality in the print.

Here is the same graphic, but showing some of the raster data, along with some wireframes.

Answer 2

I know I'm late to the party. But here is my suggestion.

1) image size
If your 550GB originals are uncompressed you should convert them to jpeg compressed tiff files. Keep them indivually (not merged). You can compress using arcgis, gdal, whatever you like. Compression will get you to around 23GB. Do not create pyramids/overviews just yet. To compress you can use any gis program you like, but I like using gdal so the command is basically this:

gdalwarp -multi -wm 512 --config  GDAL_CACHEMAX 512 -co compress=jpeg -co tiled=yes -co jpeg_quality=70 -co PHOTOMETRIC=YCBCR -co INTERLEAVE=PIXEL uncompressed.tif compressed.tif

You can easily make a bat file that goes through all your uncompressed tiffs. I like to use gdalwarp to compress my imagery instead of the usual gdal_translate, because it is faster (using multi option for multi-core, and -wm for plenty of memory).

2) handling as a single image
You can create a "virtual" mosaic using gdal vrt format. This is compatible with arcgis, qgis, mapserver, etc. Not sure about global mapper and mapinfo. The .vrt format is just a single xml file listing your images. That's a single command to create:

gdalbuildvrt nameofmosaic.vrt compressed_tif_folder\*.tif

This file is a few kb in size.

3) speeding visualization
You have to build pyramids/overviews. Just use your prefered software for this. Keeping with gdal tools you can:

gdaladdo -ro -r average --config COMPRESS_OVERVIEW JPEG --config JPEG_QUALITY_OVERVIEW 70 nameofmosaic.vrt 2 4 8 16 32 64 128

This will take a loong time. Be ready to wait 2 to 3 days of non-stop processing.

4) using the mosaic
Load the virtual mosaic in your gis program. It will be fast because it's reading the overviews which are just in 1 file like an ecw. When you zoom in to the real resolution of your images, then only the few visibile ones from the compressed images will be read, and that is really fast too.

5) handling your no data areas that show black
You have 3 solutions for this: i) use a file format that handles nodata, which is going to be complicated; or ii) use an alpha band or iii) a mask file. You can create an alpha band automatically on step 2 by telling GDAL you want nodata areas to be in the alpha band - you just add the option -addalpha:

gdalbuildvrt -addalpha nameofmosaic.vrt compressed_tif_folder\*.tif

The issue with alpha bands is that they compress badly. So your overviews are going to be larger. If that's ok with you then you're done.

If you want to create a mask file, then it is a bit more complicated. And I find this does not fit in the present question.

So hope this helps. You can find info on gdal tools by googling. Lots of interesting stuff around.

Answer 3

550gb of input TIF data is easily handled by a single ECW file. We have many customers compressing much larger datasets than this so please do not think the format is not capable in this area.

Your strategy of splitting the project into small tiles to minimize null area is also a good approach to take with the current format version as it will reduce the compression time

Your example includes reference to unsigned 16bit input data. I would recommend rescaling to 8bit if possible (depending on your requirements)

Please elaborate on why you were unable to process your project using IMAGINE or ERMapper as without this information I cannot help you. Or better yet contact the local support team

Be advised that by using the ESRI Mosaic Dataset format, what the answers above do not mention is the requirement to generate the pyramid/overview layer. Without it performance will suffer considerably. Its likely you could create the ECW equivalent files in the same amount of time but have improved image quality and significantly smaller output storage requirements.

Answer 4

Although its clearly better to use one of the other options mentioned you could try the following:

gdalbuildvrt index.vrt *.tif
gdal_translate -of "GTiff" -co "COMPRESS=LZW" -co "TILED=YES" -co "BIGTIFF=YES" index.vrt out.tif

This builds a GDAL virtual format and then convert to a single GeoTiff.

Answer 5

That sounds pretty familiar to me, we also produce large single ECW Files out of 500 too 1TB of TIF Files. But I wouldnt last on ArcGIS (ArcObjects and the Geoprocessing Engine), since it is not able to mosaic this amount in a reliable way. If you want to stay in the ESRI World I would recommend to mosaic chunks of about 50 GB or even smaller at once to a Raster Dataset stored in a File Geodatabase. The mosaic Tool tends to crash after a while, so its a good idea to let ArcGIS free memory after some GigaBytes.

Another possibility is to Use an Enterprise or Workgroup SDE Geodatabase. With SDE you get the old fashioned SDE Command Line Tools, which are built up on a robust C++ architecture other than the unreliable ArcObjects stuff. With the "sderaster -o mosaic..." command, you can mosaic to a RasterDataset until your Database store is full. There are also Commands to build up pyramids an statistics for the RasterDataset, otherwise it isnt very useful, because most Clients cannot hold the imagery in memory when reading it, as blah238 mentioned above. But Pyramids (in fact spatial indexing) should solve this problem.

But these solutions do not help you with MapInfo for sure. You mentioned you have tried ERDAS Mapper already. Thats also the Tool I would prefer. We already mosaicked 16000 TIF files, each 50 MB in size together, which are 800GB.Then we compressed it to a single ECW with a compression ratio of 1:20 which resulted in a 30GB ECW File. I am wondering that this is not working for you...

At least the whole process was running on a single core Pentium 4 1,6 GHz with 2 GB RAM, so the hardware should not be the problem. We are using Windows Server 2003 (or another server Operating System) because it better uses the harware ressources. Keep in mind, that the whole compression process needs a lot of time. Our machine was working about 5 weeks on that single file, and because it sometimes also crashed we had to do this several times, but in the end we got our ECW File.

I am not knowing another system or mechanism to store big amounts of rasters in a mostly vendor neutral way. All above mentioned ways are very ESRI specific. At least with Oracle RASTER and a pretty similar implementation in PostGIS, there are two databased Variants which are, well also not vendor neutral but opened through the SQL/MM Interface.