[GIS] Issue with Mosaic to New Raster tool

arcpygeoprocessingraster

I am having a major issue with the Mosaic to New Raster tool when importing ASC files. It processes the first few files but not all of them. I thought it was a space issue and tried to split the input but this doesn’t help.

For ex the output when I run the tool directly and in the python window is out1

enter image description here
[code]
Direct tool

arcpy.MosaicToNewRaster_management("L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8120000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8121000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8122000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8123000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8124000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8125000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8126000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8127000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8128000_1k_1m_DEM_ESRI.asc;L:/Data_Admin/temp_dataset_storage/LiDAR/FinalSupply/Atherton_Tolga_Walkamin_Mareeba_Bilwon/ASCII_DEM/asc/331000_8129000_1k_1m_DEM_ESRI.asc","P:/2011/Job_106_TRC_FireHazard_Mapping/Working/data/1_Mosaic DEM","Mosaic2","PROJCS['GDA_1994_Transverse_Mercator',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]]","8_BIT_UNSIGNED","#","1","LAST","FIRST")

[/code]

Python Code – Direct in Arcmap python window

[code]

>>> List

['L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8115000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8116000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8117000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8118000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8119000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8120000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8121000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8122000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8123000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8124000_1k_1m_DEM_ESRI.asc', 'L:\\Data_Admin\\temp_dataset_storage\\LiDAR\\FinalSupply\\Atherton_Tolga_Walkamin_Mareeba_Bilwon\\ASCII_DEM\\asc\\331000_8125000_1k_1m_DEM_ESRI.asc']

>>> arcpy.MosaicToNewRaster_management(List, v1_Mosaic_DEM, 'Name',"PROJCS['GDA_1994_Transverse_Mercator',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]]", "32_BIT_FLOAT", "", "1", "LAST", "FIRST")

<Result 'P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Working\\data\\1_Mosaic DEM\\Name'>

[/code]

So both processes are loading all the required asc files into the tool but they are not being merged and weirdly the result DEM’s are in different zones of the band.

I tried to break the files up in to bands as I thought it was a space issue but actually I need each mosaic to be made of about 100 tiles. When I tried the process by splitting into bands, it just processed the first few tiles of each region as in the graphic above.

The full script I need to run is attached and produced the output below (this prompted me to try to split into bands). I have over 250GB of free space on my C: and terabytes on Z:\ (the scratch disc)

enter image description here
Best,

# Import arcpy module
import arcpy, csv
from arcpy import env

env.scratchWorkspace='Z:\\scratch'

#model_1_2.py

def model_1_2 (RN, RootPath, List, FileList):
        # Local variables:

        for File in FileList:
                List.append(RootPath+"\\"+File[0]+"_1k_1m_DEM_ESRI.asc")

        #v1_Mosaic_DEM = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\1_Mosaic DEM"
        v1_Mosaic_DEM = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Working\\data\\1_Mosaic DEM"
        MosaicDEM = v1_Mosaic_DEM+"\\"+RN
        #MosaicDEM = RN
        slpe_ = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\3_Slope of Mosaic\\slpe_"+RN
        asp_ = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\2_Aspect of Mosaic\\asp_"+RN

        #print List
        #print "\n"+v1_Mosaic_DEM+", "+MosaicDEM+", "+slpe_+", "+asp_

        ## STEP 1

        # Check out any necessary licenses
        arcpy.CheckOutExtension("spatial")

        # Process: Mosaic To New Raster (1.1)
        print "Mosaic Processing..."
        arcpy.MosaicToNewRaster_management(List, v1_Mosaic_DEM, RN, "PROJCS['GDA_1994_Transverse_Mercator',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]]", "32_BIT_FLOAT", "", "1", "LAST", "FIRST")
        print "Mosaic Created "+ str(MosaicDEM)+"\n"

        # Process: Slope (1.2)
        print "Slope Processing..."
        #arcpy.gp.Slope_sa(MosaicDEM, slpe_, "PERCENT_RISE", "1")
        print "Slope Created "+ str(slpe_)+"\n"

        # Process: Aspect (1.3)
        print "Aspect Processing..."
        #arcpy.gp.Aspect_sa(MosaicDEM, asp_)
        print "Aspect Created "+ str(asp_)+"\n"

        # Check in any necessary licenses
        arcpy.CheckInExtension("spatial")

        ## STEP 2

        # Check out any necessary licenses
        arcpy.CheckOutExtension("3D")

        # Process: Reclassify (2.1)
        print "Aspect Reclass Processing..."
        asp__r = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\4_Reclassification_aspect_slope\\asp_"+RN+"_r"
        #arcpy.Reclassify_3d(asp_, "Value", "-1 90 1;90 180 2;180 270 3;270 315 4;315 360 5", asp__r, "DATA")
        print "Aspect Reclass Created "+ str(asp__r)+"\n"

        # Process: Raster to Polygon (2.1.1)
        print "Raster to Polygon...Aspect"
        asp__v = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\atwm.gdb\\asp_"+RN+"_v"
        #arcpy.RasterToPolygon_conversion(asp__r, asp__v, "SIMPLIFY", "VALUE")
        print "Aspect Polygon Created "+ str(asp__v)+"\n"

        # Process: Dissolve (2.1.2)
        print "Disolve Processing..."
        asp__d = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\atwm.gdb\\asp_"+RN+"_d"
        #arcpy.Dissolve_management(asp__v, asp__d, "", "", "MULTI_PART", "DISSOLVE_LINES")
        print "Aspect Dissolve Created "+ str(asp__d)+"\n"


        # Process: Reclassify (2.2)
        print "Slope Reclass Processing..."
        slpe__r = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\4_Reclassification_aspect_slope\\slpe_"+RN+"_r"
        #arcpy.Reclassify_3d(slpe_, "Value", "0.00053947966080158949 5 1;5 10 2;10 20 3;20 30 4;30 1488.7900390625 5", slpe__r, "DATA")
        print "Slope Reclassify Created "+ str(slpe__r)+"\n"

        # Process: Raster to Polygon (2.2.1)
        print "Raster to Polygon...Slope"
        slpe__v = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\atwm.gdb\\slpe_"+RN+"_v"
        #arcpy.RasterToPolygon_conversion(slpe__r, slpe__v, "SIMPLIFY", "VALUE")
        print "Slope Polygon Created "+ str(slpe__v)+"\n"

        # Process: Dissolve (2.2.2)
        print "Disolve Processing..."
        slpe__d = "P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Output\\Mosaic_FINAL_deliverable\\LIDAR\\ATWM\\atwm.gdb\\slpe_"+RN+"_d"
        #arcpy.Dissolve_management(slpe__v, slpe__d, "", "", "MULTI_PART", "DISSOLVE_LINES")
        print "Slope Dissolve Created "+ str(slpe__d)+"\n"


        # Check In any necessary licenses
        arcpy.CheckInExtension("3D")


RootName = ['BiKh', 'Mba', 'Wlkm', 'TlgTnr','Ath']
for RN in RootName:
        RootPath = r'L:\Data_Admin\temp_dataset_storage\LiDAR\FinalSupply\Atherton_Tolga_Walkamin_Mareeba_Bilwon\ASCII_DEM\asc'
        List = []
        FileList = csv.reader(open('P:\\2011\\Job_106_TRC_FireHazard_Mapping\\Working\\data\\'+RN+'-AOI.csv'))
        model_1_2 (RN, RootPath, List, FileList)
        print "\n \n Completed model for: " + RN + "\n \n"

Best Answer

There is at least one bug associated with that tool. See NIM069864 if some of the input images don't have a defined spatial reference and NIM013364 if your output is a multiband raster2. In fact the whole bug list can be found on the forum website using this [direct link] which was found by searching "Mosaic to new raster" as the search and filtered by "bugs".

Related Solutions

[GIS] Inserting DATAFRAME_ELEMENT value into dynamically generated shapefile containing its extent

I'm a bit lost here, so my answer may be off base. When you say you want the ExtentID, that appears to be the name of the Dataframe you're exporting the extent of. If this is the case (which is what your code is doing to my best guess), you don't need the [x] index when printing it. :

 for el in arcpy.mapping.ListLayoutElements(mxd2, "DATAFRAME_ELEMENT"):
     ExtentIDStore = str(el.name)
     #print x, y, z <--- You don't need to print x y and z. They are unrelated at this point.
     print ExtentIDStore #used to see which item to add to the shape file

     print ExtentIDStore[x] # <- This will just print an indexed character in the string

By the time you reach the loop above, your code is only concerned with the Dataframe Element in question, so there is no need to do any indexing for storage. ExtendIDStore should be what you add to the shapefile row.

Your error which results at this line:

row.strField = str(x)

Is probably because you are not referencing a field that exists. What comes after the row reference is the name of the Field. In this case, it should be:

row.ExtentID = ExtentIDStore

File name you can get via the filePath property of the MapDocument class and the os.path.basename method, i.e.:

filepath = mxd2.filePath # Gets full path of MXD used in your dataframe loop
filename = os.path.basename(filepath) # Extracts the file name and extension
row.FileRef = filename

If you don't want the extension, use filename.strip(".mxd") to remove it.

Scale would be:

scale = el.scale # Scale of the dataframe referenced in the loop
row.Scale = scale

And that should fill in the rows you're looking for.

For the new problem, I think you are misusing the insert cursor. The second argument should be a spatial reference, which cannot be gleaned via an MXD. I would recommend creating the cursor with only one argument (the output file or dfStore). You also forgot the last (and most important) line:

rows.insertRow(row)

Without that your data never gets put into the table. Put that down at the end of the row loop so that the data is actually inserted after all the variables are set.

For more information about accessing Dataframe, layer, and map properties (scale, name, current position), check out the mapping module. The scale property of the Dataframe class can be seen on the classes page of the same name.

[GIS] Mosaic to New Raster tool fail

This is a common ArcGIS problem. I believe as you zoom-in some of the areas that seem not to be there begin appearing as if it was a rendering issue...

I think, although it has never been confirmed, that it most often happens when mosaicing tiles from File Geodatabase (FGDB) sources and outputting the mosaic raster to the same or other File Godatabase.

My crude workaround is to import the individual source tiles into a File Geodatabase but when mosaicing output the result into a folder as ESRI Grid. When complete import back into FGDB.

I have most success using this approach but also had it fail a few times. When it does, I try mosacing from individual ESRI Grid files and output to an ESRI Grid. When all fails i try the reverse, and so on.

Sounds strange but works for me.

Best Answer

Related Solutions

[GIS] Inserting DATAFRAME_ELEMENT value into dynamically generated shapefile containing its extent

[GIS] Mosaic to New Raster tool fail

Related Question