Clipping Rasters from Mosaic Dataset in Original Coordinate System?

Question

I am building a large mosaic dataset that contains several hundred rasters, which are spread across three UTM zones. My goal is to be build a tool in ArcMap (most likely within model builder) that will clip the mosaic dataset by a user-selected shapefile, which would then be exported for use in another software program.

The secondary software utilizes 8 digit UTM coordinate systems (i.e., has the UTM zone number in front of the easting). As such, I've made custom coordinate systems covering the study area which have the UTM zone number included in front of the false easting (e.g., UTM zone 14N would have a false easting of 14500000m and so on).

The rasters in the mosaic dataset are already projected in these custom coordinate systems, but the mosaic dataset itself is in WGS84 World Mercator, although that could be changed if it would help.

Whenever I clip the mosaic dataset via feature class, the output raster is in the coordinate system of the mosaic dataset (WGS84 World Mercator), while I would like it to be in the orginal coordinate system of the source raster.

I'm looking for a way to either retain the coordinate system of the source raster(s) when clipping or automate a tool to project them in the correct coordinate system based on their location or UTM zone. I have added the correct UTM zone of each source raster in the mosaic dataset's attribute table, but don't know where to go from here.

I'm not very familiar with Python or programming in general, so I would prefer to create a tool in model builder, but I am not opposed to making a python script if that's the only potential method.

One idea, although I'm not sure how to execute it, would be to run an If Then type statement on the location of the feature class used for clipping. I.e., if the feature class is within a shapefile representing Zone 15, then use 'X' projection. If the clipping feature class is within a shapefile representing Zone 14, then use 'Y' projection, etc.

After @Michael's comments, I'll work on integrating that script. Looking at this link http://resources.arcgis.com/en/help/main/10.1/index.html#//002w00000022000000, I am eventually trying to build a similar true/false type logic. The goal would be to run a script testing which UTM zone the clipping shapefile is in, and then using the results of that statement to set the coordinate system for the extraction process.

Any thoughts on how to accomplish that? I was thinking testing whether the extraction shapefile's centroid is contained within a shapefile representing an entire UTM zone. IF I had 3 shapefiles (one for each UTM zone in the study area), the extraction shapefile could be tested somehow to see in which shapefile it is contained, thus providing a "true" response in that instance. Something along the lines of "if the extraction shapefile's centroid is within the zone 14 shapefile, use coordinate system A. If the extraction shapefile's centroid is within the zone 15 shapefile, use coordinate system B" etc.

Answer 1

Set the Output Coordinate System environment setting when clipping the raster using Extract by Mask with an input shape file. If you are implementing this in python you should be able to grab the coordinate system from the shape file using the describe method.

Projecting/converting WGS84 World Mercator to WGS84 UTM zone 14N does not require a parameter transformation as they are based on the same geographic coordinate system, if they were dissimilar then this would not work.

Here is a script that will clip a raster mosaic to a shapefile and match the output coordinate system to the shapefile:

import arcpy, sys

InRasterMosaic = sys.argv[1]
InShapeFile = sys.argv[2]
OutRaster = sys.argv[3]

if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.AddMessage("Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")
else:
    arcpy.AddMessage("Unable to get spatial analyst extension")
    sys.exit(0)


desc = arcpy.Describe(InShapeFile)
SR = desc.spatialReference
arcpy.env.outputCoordinateSystem = SR

ExRas = arcpy.sa.ExtractByMask(InRasterMosaic,InShapeFile)
ExRas.save(OutRaster)
arcpy.CheckInExtension("Spatial")

I have tested this and it works! Have a read of Adding a Script Tool, follow those instructions and you can put this into a model.

If the input is in WGS (DD) as stated:

import arcpy, sys, os

InRasterMosaic = sys.argv[1]
InShapeFile = sys.argv[2]
OutRaster = sys.argv[3]

if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.AddMessage("Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")
else:
    arcpy.AddMessage("Unable to get spatial analyst extension")
    sys.exit(0)

# from the geographic extent, get the centre
# of the extent, only the X is required
desc = arcpy.Describe(InShapeFile)
Ext = desc.Extent
CentX = ( Ext.XMin + Ext.XMax ) / 2

# I'm not familiar with North American projections
# so you'll have to put the numbers in for each zone
if CentX < Zone14max:
    ProjSR = arcpy.SpatialReference(32614)
elif CentX < Zone15max:
    ProjSR = arcpy.SpatialReference(32615)
# continue for the required spatial references
else:
    arcpy.AddMessage("Unable to establish coordinate system")
    arcpy.CheckInExtension("Spatial")
    sys.exit(-1)

arcpy.env.outputCoordinateSystem = ProjSR

ExRas = arcpy.sa.ExtractByMask(InRasterMosaic,InShapeFile)
ExRas.save(OutRaster)

arcpy.CheckInExtension("Spatial")

This will need you to make some edits, but will basically do what you're after. Using the centroid of the extent work out which zone it's likely to be in and then use that for the output coordinate system. The numbers are the SRID of the output coordinate system like this one http://spatialreference.org/ref/epsg/32614/