[Gdal-dev] Need a sanity check on processing imagery

[Stephen Woodbridge]

Hi all,

I'm building a process to extract 3.75 x 3.75 min GTif tile images from 
HUGE CCM MrSID files and I am hoping I can prevail on you to review this 
and let me know if I'm doing some horridly wrong.

The MrSID files are in UTM and we want to project all the tiles into 
geographic (latlong) projection.

$ gdalinfo -nomd naip_1-1_1n_s_ut057_2006_1.sid
Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Size is 95310, 41990
Coordinate System is `'
Origin = (373650.000000,4588120.000000)
Pixel Size = (1.00000000,-1.00000000)
Corner Coordinates:
Upper Left  (  373650.000, 4588120.000)
Lower Left  (  373650.000, 4546130.000)
Upper Right (  468960.000, 4588120.000)
Lower Right (  468960.000, 4546130.000)
Center      (  421305.000, 4567125.000)
Band 1 Block=1024x128 Type=Byte, ColorInterp=Red
   Overviews: 47655x20995, 23828x10498, 11914x5249, 5957x2625, 
2979x1313, 1490x657, 745x329, 373x165, 187x83, 94x42, 47x21, 24x11, 
12x6, 6x3
Band 2 Block=1024x128 Type=Byte, ColorInterp=Green
   Overviews: 47655x20995, 23828x10498, 11914x5249, 5957x2625, 
2979x1313, 1490x657, 745x329, 373x165, 187x83, 94x42, 47x21, 24x11, 
12x6, 6x3
Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue
   Overviews: 47655x20995, 23828x10498, 11914x5249, 5957x2625, 
2979x1313, 1490x657, 745x329, 373x165, 187x83, 94x42, 47x21, 24x11, 
12x6, 6x3

I'm getting the utm zone from here and it just dawned on me that I can 
do something like:  -s_srs naip_1-1_1n_s_ut057_2006_1.prj instead of 
reading the file myself. Right?

$ cat naip_1-1_1n_s_ut057_2006_1.prj
PROJCS["NAD_1983_UTM_Zone_12N",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
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",0.0],
PARAMETER["Central_Meridian",-111.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]

I have a perl script that inverse projects the corners of the sid back 
to lat, lon then iterates over the file generating commands like this 
for each to the 3.75x3.75 min tiles:

$ gdalwarp -s_srs +init=epsg:26912 -t_srs EPSG:4326 -te -111.9375 
41.4375 -111.875 41.5 -rc -wm 250  -co "TILED=YES" 
naip_1-1_1n_s_ut057_2006_1.sid B-111_41_129.tif

Some of the tiles will have to be mosaiced into a common tile from 
multiple CCM files along their edges. Some of the options above imply 
that I want a new file and NOT a mosaic. Would that be the -t_srs 
EPSG:4326 and -co "TILED=YES" options?

$ gdalinfo B-111_41_129.tif
Driver: GTiff/GeoTIFF
Size is 5420, 5420
Coordinate System is:
GEOGCS["WGS 84",
     DATUM["WGS_1984",
         SPHEROID["WGS 84",6378137,298.2572235629972,
             AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0],
     UNIT["degree",0.0174532925199433],
     AUTHORITY["EPSG","4326"]]
Origin = (-111.937500,41.500000)
Pixel Size = (0.00001153,-0.00001153)
Metadata:
   AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (-111.9375000,  41.5000000) (111d56'15.00"W, 41d30'0.00"N)
Lower Left  (-111.9375000,  41.4375034) (111d56'15.00"W, 41d26'15.01"N)
Upper Right (-111.8750034,  41.5000000) (111d52'30.01"W, 41d30'0.00"N)
Lower Right (-111.8750034,  41.4375034) (111d52'30.01"W, 41d26'15.01"N)
Center      (-111.9062517,  41.4687517) (111d54'22.51"W, 41d28'7.51"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue

Does this look like a reasonable process? Or is it flawed?
Can I do anything to speed this up? make -wm larger?
Can I extract multiple file tiles without having to open the file each 
time? Can that be done in Perl with the Geo::GDAL module?

I would appreciate and help and advice before a loose this monster on a 
few TB of files.

Many thanks in advance,
   -Steve