[Gdal-dev] I think I need help with warping.

[Frank Warmerdam]

Kyler Laird wrote:
> I'm trying to combine lots of map images.  Here are a couple of small edge
> samples with a lot of overlap.
> 	http://aviationtoolbox.org/tmp/translate1.tif
> 	http://aviationtoolbox.org/tmp/translate2.tif
> 
> I'm confuesd by the need to rotate/warp these in order to make them line up to
> each other.  They appear to be in the same projection (Lambert Conformal Conic
> 2SP) but the latitude_of_origin, central_meridian, and pixel size values are
> slightly different.  I'm still working on understanding the meaning of those
> differences.

Kyler,

If they have different parameters for the projection then they are different
coordinate system with the same underlying projection method.  You need to pick
a coordinate system and resolution you want to operate on country wide (or perhaps
split the country into a few sections each in their own coordinate system.

> However...while I'm pursuing enlightenment, it would help me immensely if
> someone would make some suggestions on combining these images.  Eventually I
> want to go across the entire ConUS and I need to redo it every 28 days so I
> want to do this well - not just manually slop together something that looks
> about right.

The trickiest part of this is trimming off all the stuff that *isn't* part of the
image itself.  I have looked at the Atlanta north and south map chunks, and it
doesn't look trivial to me to cut away all the surround information.  The map data
area's don't seem to fall on clean cutlines in lat/long or the projection that is
used.

I have made a crack at doing this and this is the process I used:

  o In order to trim off the map surround, I collected a polygon cutline for
    the area *off the map* in OpenEV, and saved it as a shapefile.

  o I wrote a script in python, but using OpenEV as well as GDAL instrisics
    to trim off the surround area. This was accomplishing using the OpenEV
    polygon rasterization intrinsics.    The script is attached (clip_raster.py)
    The clip_raster.py needs to be able to update the file in place, so I made
    a copy of the source file using gdal_translate first ... the copy is by
    default uncompressed and should be updatable in place with GDAL.  The
    original might have been ok too, but I didn't want to mess with it.

    The result of running clip_raster.py is that the areas under the polygon
    (presumably the non-map areas) are masked out with the value 255 which seems
    to be used for some no-data edging in the original tiff files.

    % gdal_translate Atlanta71North.tif at_north_wrk.tif
    % clip_raster.py at_north_wrk.tif mask_at_north.shp

  o Next, I converted the files to RGB so they can be merged effectively.  I
    just used the stock pct2rgb.py script in gdal/pymod for this.

    % ~/gdal/pymod/pct2rgb.py at_north_wrk.tif at_north_rgb.tif

  o Repeat above sequence for other files, in my case the atlanta south file.

  o Now I needed to decide on an output projection for the mosaic file.
    For this effort I just used the projection the atlanta files were
    already in.  Copied from the gdalinfo output it looks like
    this (which I saved to a file 'atlanta.wkt'):

PROJCS["unnamed",
     GEOGCS["NAD83",
         DATUM["North_American_Datum_1983",
             SPHEROID["GRS 1980",6378137,298.2572221010042,
                 AUTHORITY["EPSG","7019"]],
             AUTHORITY["EPSG","6269"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4269"]],
     PROJECTION["Lambert_Conformal_Conic_2SP"],
     PARAMETER["standard_parallel_1",38.66666666666667],
     PARAMETER["standard_parallel_2",33.33333333333336],
     PARAMETER["latitude_of_origin",34.10833333333333],
     PARAMETER["central_meridian",-84.50000000000003],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

  o Next, I needed to create a mosaic file with extents in the selected
    output projection suitable to hold all the inputs.  To compute the
    appropriate bounds and create the empty mosaic file I created a script
    (attached as create_mosaic.py).

    % create_mosaic.py ./atlanta.wkt mosaic.tif at_north_rgb.tif at_south_rgb.tif

  o The mosaic file now exists, but is empty. I now use gdalwarp to "copy" each
    file in.  At this step it is necessary to know what color 255 translated into
    in the pct2rgb process.  It seems the value mapped to different colors in the
    two files I got.  I used OpenEV to inspect the output image from the pct2rgb
    process, but you should also be able to inspect the last entry in the input
    file color table.  I got the rgb values 206,41,33 for the north file, and
    8,107,8 for the south file.

    % gdalwarp -srcnodata '206 41 33' at_north_rgb.tif mosaic.tif
    % gdalwarp -srcnodata '8 107 8' at_south_rgb.tif mosaic.tif

This produced a pretty decent mosaic.  Some issues are:

  o You are unlikely to be able to TIFF format for the mosaic since even in this
    small case of two parts of a map, the resulting file was 334MB and the limit
    for TIFF is 4GB.

  o The quality of the match between mapsheets wasn't great.  I would guess the
    georeferencing of the individual files can be out by as much as 240 meters or
    so.  Still ... no bad.

  o For some maps the edges are quite curved, meaning the polygon cutline needs to
    have alot of vertices along the curve.

  o The cutlines and "nodata values" may not be reusable when the maps are
    regenerated ... but it is hard to be sure since I don't know the process
    for generating them.

A proof of concept, though not a problem free process. Some of the steps could
be boiled into a more automated script.

Best regards,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

-------------- next part --------------
#!/usr/bin/env python

import gdal
import osr
import gview
import gvalg
import sys

# =============================================================================
def Usage():
    print 'Usage: clip_raster.py <raster> <mask_in>.shp'
    sys.exit( 1 )

# =============================================================================
#
# Program mainline.
#

if len(sys.argv) != 3:
    Usage()

raster_filename = sys.argv[1]
mask_filename = sys.argv[2]

# =============================================================================
# Instantiate a GvRaster for the file in update mode.

raster_fd = gdal.Open( raster_filename, gdal.GA_Update )
raster = gview.GvRaster( dataset = raster_fd )

# =============================================================================
# Load the vector layer. 

mask = gview.GvShapes( shapefilename = mask_filename )

# =============================================================================
# Apply masking.

gvalg.rasterize_shapes( raster, mask, 255 )

-------------- next part --------------
#!/usr/bin/env python

import gdal
import osr
import sys
import osr

# =============================================================================
def Usage():
    print 'Usage: create_mosaic.py <coordinate_sys> <mosaic> <raster> ...'
    sys.exit( 1 )

# =============================================================================
# Program mainline.
#

if len(sys.argv) < 4:
    Usage()

srs_arg = sys.argv[1]

mosaic = sys.argv[2]
files = sys.argv[3:]

# =============================================================================
# Instantiate the coordinate system.

target_srs = osr.SpatialReference()
if target_srs.SetFromUserInput( srs_arg ) != 0:
    print 'Failed to create SRS.'
    sys.exit( 1 )

print 'Generating Mosaic in Coordinate System:'
print target_srs.ExportToPrettyWkt()

# =============================================================================
# Process all files.

targ_xmin = None
targ_ymin = None
targ_xmax = None
targ_ymax = None

for filename in files:

    fd = gdal.Open( filename )

    # Read projection, and create transformer to target SRS. 
    file_srs_wkt = fd.GetProjection()
    file_srs = osr.SpatialReference()
    file_srs.ImportFromWkt( file_srs_wkt )
    
    ct = osr.CoordinateTransformation( file_srs, target_srs )
    
    # Setup bounds in source srs.
    gt = fd.GetGeoTransform()
    
    ulx = gt[0]
    uly = gt[3]
    
    lrx = gt[0] + fd.RasterXSize * gt[1] + fd.RasterYSize * gt[2]
    lry = gt[3] + fd.RasterXSize * gt[4] + fd.RasterYSize * gt[5]
    
    corners = [ (ulx, uly), (ulx, lry), (lrx, uly), (lrx, lry) ]
    
    for corner in corners:
        (x,y,z) = ct.TransformPoint( corner[0], corner[1] )
        if targ_xmin is None:
            targ_xmin = x
            targ_xmax = x
            targ_ymin = y
            targ_ymax = y
        else:
            targ_xmin = min(targ_xmin,x)
            targ_xmax = max(targ_xmax,x)
            targ_ymin = min(targ_ymin,y)
            targ_ymax = max(targ_ymax,y)

if targ_xmin is None:
    print 'Failed to get extents, stopping.'
    sys.exit( 1 )

# =============================================================================
# What is the target file size and geotransform?
#
# For now we will hardcode 60m x 60m pixels.

x_pixel_size = 60.0
y_pixel_size = 60.0

xsize = (targ_xmax - targ_xmin) / x_pixel_size + 1
ysize = (targ_ymax - targ_ymin) / y_pixel_size + 1

targ_gt = (targ_xmin, x_pixel_size, 0.0,
           targ_ymax, 0.0, - y_pixel_size)

# =============================================================================
# Create target RGB "mosaic" file.  

driver = gdal.GetDriverByName( 'GTiff' )

fd = driver.Create( mosaic, xsize, ysize, 3, gdal.GDT_Byte, [ 'TILED=YES' ] )

fd.SetProjection( target_srs.ExportToWkt() )
fd.SetGeoTransform( targ_gt )