Saturday 24 November 2018

arcgis desktop - How to workaround large mosaic process which is failing?


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.


enter image description here 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---


enter image description here


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



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.


enter image description here


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


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...