Scripting to load only selected rasters based on an index shapefile

I realize I still haven’t gotten around to Part 2 of that other sequence yet. I will, but I hit a little nuisance over the weekend at work…well, it’s always been a nuisance, but a minor one. As we move to higher resolution imagery, pixel counts rise, file sizes or numbers increase, and I foresee that this is a nuisance that will grow, not go away, so I thought I’d see if I could make my life easier and my work more efficient–my time better spent, my employers’ money better invested–if I could figure out a way to load only the files I need to cover a certain area of interest (market area, administrative boundary,…?) without manually selecting them from a folder containing 6000 files. For whole rows, it’s not difficult…just select a range, but for some cases that’s inconvenient, like if you want a slightly irregular chunk down the eastern edge cut into six .ecws, stair-stepped along the eastern corner of the county…then you want sequences of uneven length from the right side of each row, also of uneven length, and you have to manually select these out of your window listing all the 6000 .tifs for the county to load into Global Mapper because reprojecting the whole county would tie up your machine for a length of time that you might like to reduce…

Where was I? Oh yes. It’s pretty easy to pop into ArcMap, overlay the two index grids, and Select By Location to figure out which .tifs I really need. I know from noodling around that it’s pretty easy to create a Search Cursor in Python to read all the filenames out of a field in my shapefile and spit them into a text file. This was a  n early babystep–I could at least then copy and paste 20 or so at a time from the text file into the Open>File dialog box, without having to manually look and remember or write down “10682 to 10473, 9986 to 9840, …” for each row. Just blithely cut and paste, 20 at a time. For 6 sids, at what, 300, 400 .tiffs at a time, maybe more? Yeah. Better, but not great.

So then the next step was figuring out how to automate the ingestion into Global Mapper of this information. Ideally, if ArcMap could talk to Global Mapper, I could have ArcMap tell Global Mapper which tiles I want to load.

Incidentally, I found that having Global Mapper import files, reproject, and export via script it works better and faster on my machines at work than it does if I have to wait for the display to refresh during the reprojection. And that was a direct offshoot of this.

I had noticed in Global Mapper’s File drop-down an option to “Run Script”, so I thought, hmm. If it understands Python, that would be handy, since I’ve already played with that a bit. It doesn’t. I looked into Global Mapper’s forum and Scripting Language Reference, and discovered in the process that the workspace (.gmw) is actually just a script that loads the files and sets the projection and display characteristics. It’s really pretty simple when you stare at it for a while. And exporting to a raster format is really pretty much just as easy as importing and setting up the workspace, except it says “export” instead of “import”.

So why not set up a tool in ArcMap to export a Global Mapper script (.gms), I thought. It appears the scripting in Global Mapper is such that I could probably set the whole thing up to run directly at a command from a tool integrated into ArcMap, but for now I figured I’d just run the scripts from IDLE and code paths directly in a text editor or shell, rather than deal with GUI inputs.

Anyhoo, so there you have it. I take two shapefiles and overlay them to find which features in the one layer intersect a certain feature or set of features from the other layer. Then I figured out what the Global Mapper script would need to say, which was basically:

GLOBAL_MAPPER_SCRIPT VERSION=1.00

followed by

IMPORT “C:\GIS\Blog\GeoTIFFs\” + file_name over and over.

This is what my Python script would write. Then if I wanted to, I could project it, say, from State Plane 6 to UTM 11 N, and export it from the script by adding a couple lines saying so after the import:

IMPORT FILENAME=”AOI.shp”

DEFINE_PROJ PROJ_NAME=”UTM_ZONE11_NAD83″
Projection UTM
Datum D_NORTH_AMERICAN_1983
Zunits NO
Units METERS
Zone 11
Xshift 0.000000
Yshift 0.000000
Parameters
END_DEFINE_PROJ

LOAD_PROJECTION PROJ_NAME=”UTM_ZONE11_NAD83″

EXPORT_RASTER FILENAME=”AOI_2011_3in.ecw” TYPE=ECW LAYER_BOUNDS=”AOI.shp” GEN_PRJ_FILE=YES TARGET_COMPRESSION=”2″

So, for simplicity and speed on my part, I would write a Python script that would take a shapefile as an input, create the header and import statements for the tiffs. To get this shapefile, I would do the selection in ArcMap and save the result. (Ultimately I could just create a layer and use the layer from a dialog box within ArcMap.) The shapefile I want it to read is cells on the tiff index grid selected because they intersect a selected grid cell on the coarser grid:

And the link between this shapefile and the Global Mapper script would be this Python script:

Then I can just open up Global Mapper, select Run Script, navigate to and select the .gms, and then hit go.

 

Leave a Reply

Your email address will not be published.