0

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.

4
  • Arbitrarily setting the coordinate system is dangerous, you cannot accommodate all the possible coordinate systems that a shape file can be in. Make the user pick one or force them to set it correctly. I know what you're trying to do but I am saying that road only leads to pain! Jun 2, 2014 at 0:26
  • I apologize for not clarifying, but the clipping shape file would be generated in Google Earth by the "user" downstream, so it would be in a known (or at least "normal" coordinate system). All I would need to do is automate a way for that shapefile to be located within one of the three UTM zones in the study area, which would then provide the basis for selecting the eventual output coordinate system.
    – Eric
    Jun 2, 2014 at 0:42
  • Oh, that's different! That would make it not so arbitrary Would it be geographic, mecator, projected? Jun 2, 2014 at 0:45
  • I just tested converting a polygon from Google Earth to ArcMap, and it was in geographic WGS84. Here are the full details: GCS_WGS_1984 WKID: 4326 Authority: EPSG Angular Unit: Degree (0.0174532925199433) Prime Meridian: Greenwich (0.0) Datum: D_WGS_1984 Spheroid: WGS_1984 Semimajor Axis: 6378137.0 Semiminor Axis: 6356752.314245179 Inverse Flattening: 298.257223563
    – Eric
    Jun 2, 2014 at 0:52

1 Answer 1

2

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/

5
  • Thank you for the reply Michael. That solution won't fulfill the automation requirement though (unless I am missing something?) Would there be a way to first project the clipping shape file into the desired 8 digit UTM coordinate system based on its location? This is assuming that the shape file would be in the "incorrect" coordinate system and would need to be converted to my custom one.
    – Eric
    Jun 1, 2014 at 23:32
  • Considering you would be grabbing the coordinate system from the shape file it would be up to the user to ensure that it's good, however you can test it and return an error if it's not one that is recognized. You simply cannot account for every possible input coordinate system that the shape file can be in. The output coordinate system environment tells the geoprocessor to produce results into that coordinate system so your extracted raster should match the shapefile. Jun 1, 2014 at 23:38
  • I edited the original question to include my response to your reply, thank you.
    – Eric
    Jun 2, 2014 at 0:20
  • Thank you very much sir, I'll work with that script and hopefully I can figure it out!
    – Eric
    Jun 2, 2014 at 1:19
  • You're welcome. You should only need to put in the values for your output coordinate systems, maximum geographic coordinate and SRID. I think zone 14 goes up to 96 West based on this page geoinfo.nmt.edu/publications/maps/gps/home.html but I'm not certain so I left it as an unset variable. Jun 2, 2014 at 1:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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