Rick -- this is a restatement of what I sketched out on the board earlier as one way of quantifying the boundary behavior:

1. Calculate the difference between the fused ASTER/SRTM layer and the CDEM layer. I'll call this set of differences D, for short. To the extent that we regard the CDEM as "topographic truth", then D basically represents the spatial pattern of fused DEM "errors". I think it's a reasonable null hypothesis that the distribution of D is invariant with respect to the 60N boundary, and we can use this as a starting point for statistically testing for the presence of boundary artifacts.

You said you already did this for a subset of the respective layers using the R raster package, but only after manually forcing the extent of one raster to match the other. As we discussed, verify that the magnitude of this artificial shift was indeed insignificant relative to the cell resolution, as you thought.

2. For exploratory purposes, perhaps compare the distribution of D values south of 60N (SRTM zone) with the distribution of values north of 60N (ASTER zone), taken for some sufficiently large set of randomly selected pixels. This won't tell you about boundary artifacts per se, but I think it would paint a useful picture of how ASTER >60N and SRTM <60N each separately compare to GDEM.

• below 60N
• above 60N
• exactly across the 60N boundary

You can then use these values to test e.g. a null hypothesis that the difference in D between pairs of N-S adjacent pixels that lie across the 60N boundary is the same as the difference between N-S adjacent pairs that do not cross the boundary. I imagine it will be good to test for differences in a few different properties of the distributions (mean? variance? other aspects of shape? something about extremes?). Beyond that, perhaps there is even a specific literature on relevant DEM fusion artifacts that bears on this? This may be exactly the kind of thing worth passing by John Gallant for feedback.

4. As we discussed, go back to look at (and report on in an appropriate Redmine issue) the details of your ASTER-SRTM merge operation. Did you simply splice the data at 60N, such that the fused layer contains unadulterated SRTM values south of the boundary, and unadulterated ASTER values to the north? Or was some sort of resampling done that has already resulting in smoothing or some such at the boundary, which probably means expanding the scope of the analysis outline above?

• DEM aggregated to 1km
• Derived terrain layers (slope, aspect, etc) at 90m
• Derived terrain layers at 1km

Everything you're doing should be fully scripted, so this shouldn't imposed any extra burden, but I just wanted to mention it.

Overlapping with Jim's notes (aboveon our discussions today:

I have created two images corresponding to southern British Columbia, bounding the 60 degree North Latitude line: the first is extracted from the ASTER/CGIAR SRTM image mosaic, the second is extracted from the Canada CDEM image mosaic.

Next, I created an R script, "filterDiffRasterImage.r", that 1) extracts a narrow (<= 100 pixels) subimage along the 60 degree Latitude line from each of these two images, using the R 'raster' package, and 2) subtracts the CDEM subimage from the ASTER/CGIAR SRTM image, creating a 'difference image'.

The next steps:

1) Resolve a possible mis-alignment between the CDEM image and the mosaic terrain image;

2) is to employ this script and the CDEM terrain map to statistically compare the difference between CDEM and each component terrain mosaic type (ASTER GDEM and CGIAR SRTM), with the difference in the mosaic image region along the boundary between image types. The object of this step is to characterize and compare the distribution of the elevation pixel differences along the image mosaic border region with the distribution of differences in the component images (ASTER, CGIAR/SRTM) in areas away from the border. Ideally, the distributions will be similar.

Jim R and I discussed a procedure for accomplishing step 2). I plan to work on this for the remainder of Tuesday and Wednesday.

Rick -- I just noticed your r-sig-geo exchange with Hijmans, and I want to expand on the associated issue I raised above about whether your raster 'shifting' approach is okay.

First, I think Robert's estimate of the error is off. For starters, if your artificially imposed shift is 0.00036 degrees (which seems to be the case based on what you reported in that post), and if the cell resolution were indeed 0.0083333 as Robert speculated, then the error would be only ~4% of the cell resolution (not 1/4). But moreover, the 3" data you are working with actually has a resolution of 0.0008333 degrees, meaning the shift is ~40% of a cell width. Note that the worst possible offset would be 50% of a cell width, because the clipped rasters won't be misaligned by more than half a pixel -- otherwise it would've clipped different pixels! So at least relative to cell resolution, this is not what I was hoping for when I said to confirm that the amount of shift is insignificant. Doesn't mean it's a fatal problem, but I don't think it can be dismissed without further probing.

As I described yesterday, and as Robert similarly indicated in his reply, the alternative is to resample the CDEM to match your fused DEM. Actually, I suppose what you did by slightly shifting the raster really equates to nearest-neighbor resampling, but now we're talking bilinear (or cubic, etc) resampling. Unfortunately, this raises its own issues because it actually alters the data values in some way. It's not completely clear to me which is the lesser of two evils, but my suggestion in cases like this is to try it both ways and see if it even matters in the first place (my guess is, probably not). As always, this should be simple enough as long as you are setting everything up in a scripted manner.

Jim:

I think that Robert pointed me in the correct direction, but I took some steps to quantify the error and test it's effect
on the difference image results by examining the CDEM and merged ASTER/SRTM images at 'close range'.

The 'shift' issue occurs because the pixels in these two terrain images have slightly different 'origins'. The difference
is less than one pixel, but a different amount, in each dimension. By superimposing the images and increasing the magnification
on the lower left corner of the image pair and using the mouse to obtain the corner coords in decimal degrees, I made this estimate:

Origin (LLC) coordinates

CDEM         ASTER/CIGAR merge    difference/offset
                                       degrees arc-sec    meters
X    -135.9882       -135.9881          .0001   0.33     ~10 meters @ 58 deg N
Y      58.4653         58.4648          .0005   1.8      ~50 meters

As you say, the offset is a bit larger than Robert estimated, esp in the Y dimension.

So, I created a resampled image in which (using the images from which I obtained the above coordinates),I sampled the
ASTER/CGIAR image into the coordinate space of the (identically-sized) CDEM image, using gdalwarp:

gdalwarp -te -135.9882 58.4653 -r bilinear asterCgiarSub.tif asterCgiarSub.tif (as I remember)

Then, I extracted the same type of 10 row 'boundary image' centered on the 60 degree N lat line, as before, and performed
image difference calculation using this image and the coincident CDEM image. I will review the actual numbers when I get
to work (after the library), but I remember that the mean of the resampled difference image was slightly larger than that
of the 'offset' image. I think that the difference was, -4 vs -3 (the offset image). The Std Dev's were very close, say,
25 vs 24. I want to reconfirm these numbers, and then create separate difference images for CDEM vs Aster and CDEM vs SRTM,
to see how they compare. This should all be done by just after lunch if not sooner. I will check in the (short) R and Shell
scripts that I created to perform this task.

My feeling is that the differences in the boundary regions will not be statistically different from the 'outlying' regions,
but I will wait for the actual numbers to see for sure.

The book is on order, and should arrive at UCSB Library in 3 to 5 working days.
It will come in handy in several parts of the Fused Terrain Layer work, but especially in constructing the derived terrain layers.

Presumably "the book" refers Terrain Analysis (2000) by Wilson and Gallant? Has the state of terrain analysis not advanced much over the past decade? I think we should be aware that the types of coverages we are using (ASTER GDEM, SRTM) did not exist when the book was written, nor did some of the currently popular analytical libraries (R raster, gdal, etc.). Just a cautionary note that the information in there might be outdated.

Yes, the book IS the Wilson and Gallant book. This book is a decade old, and new techniques no doubt have been developed
that we need to review and test as we develop workflows. So, more research is in order, even as we wait for this book
to arrive next week.

Good reasons to have a look at this book:

- It is well-reviewed
- One of the authors is a recognized authority on TA, and a member of this WG
- The table of contents contains topics relevant to this project
- The UC Library system is providing it free-of-charge :}

A few comments:

Rick Reeves wrote:

By superimposing the images and increasing the magnification
on the lower left corner of the image pair and using the mouse to obtain the corner coords in decimal degrees, I made this estimate:

I'm curious, why not obtain these values by querying them from the raster programmatically, rather than doing it manually? E.g., using gdalinfo, or extent() in the raster package, etc. Also, no need to convert the offset to meters (which of course varies by latitude in the X direction anyway). I really do think the more important thing is the offset as a proportion of cell resolution.

So, I created a resampled image in which (using the images from which I obtained the above coordinates),I sampled the
ASTER/CGIAR image into the coordinate space of the (identically-sized) CDEM image, using gdalwarp:

Important: What I had said was to do the opposite, namely resample the CDEM into the space of your fused DEM. Otherwise all of your validation is based on a resampled version on the fused DEM, not the fused DEM itself.

gdalwarp -te -135.9882 58.4653 -r bilinear asterCgiarSub.tif asterCgiarSub.tif (as I remember)

Is the exact command recorded somewhere, whether as part of a script or otherwise?

I will review the actual numbers when I get
to work (after the library), but I remember that the mean of the resampled difference image was slightly larger than that
of the 'offset' image. I think that the difference was, -4 vs -3 (the offset image). The Std Dev's were very close, say,
25 vs 24.

Good to start to see some numbers, though again you'll need to re-do things after obtaining the resampled CDEM. You should then look at the correlation coefficient between the two (i.e., between the shifted CDEM and resampled CDEM). Also, I would calculate the mean/sd/range of the difference between these two CDEM rasters, rather than comparing the stats computed separately on each of them. Lastly, to clarify my earlier comments, I think unless these two different aligned CDEM rasters are so similar that it would be obviously pointless, and assuming it is trivial to just re-run your scripted validation workflow, I think it would be useful to do the same fused layer validation calculations separately using the shifted CDEM and the resampled CDEM to see if it matters.

Jim:

As I mentioned, after I completed my initial 'screen based' review of the boundary area, I intend to perform the statistical analysis programmatically, as you say.

Important: What I had said was to do the opposite, namely resample the CDEM into the space of your fused DEM. Otherwise all of your validation is based on a resampled version on the fused DEM, not the fused DEM itself.

I initially resampled the fused DEM into the CDEM space in order to see how much the elevation values changed
during the resample operation. There was almost no change in the elevation values. This result should not change if I resampled the CDEM into the DEM space, should it?

And I intend to compute the statistics on the difference between the CDEM rasters, and present the information.

BTW, here is the R script that I used to generate the edge region and difference subimages:

require(raster)
require(rgdal)

inputFirstRaster <- raster(sFirstImageName)
inputSecondRaster <- raster(sSecondImageName)

# Create extent objects used to extract raster subimges. 
# The object will be centered along the 60 degree North latitude line,
# and have varying depths (number of rows).
# The Western Canada study area runs from -135 (west) to -100 (west) longitude,
# and 55.0 to 64.00 degrees (north) latitude. 

#eTestAreaExtent <- extent(-135.2,-100.2, 59.995,60.005) # Creates a 12 row subimage
eTestAreaExtent <- extent(-135.2,-100.2, 59.997,60.001) # Creates a 5 row subimage

# Extract a sub image corresponding to the selected extent.
# Two different alternatives:
# The extract() function returns a vector of cell values, 
# the crop() function returns a complete raster* object.

vEdgeRegionFirst <- extract(inputFirstRaster,eTestAreaExtent)
rEdgeRegionFirst <- crop(inputFirstRaster,eTestAreaExtent)
vEdgeRegionSecond <- extract(inputSecondRaster,eTestAreaExtent)
rEdgeRegionSecond <- crop(inputSecondRaster,eTestAreaExtent)

# Important: In order for the image subtraction to work, the extents
#            of the two images must be IDENTICAL. I used ArcMap GIS Raster Crop By Mask
#            to create subimages with identical extents. 

rDelta <- rEdgeRegionFirst - rEdgeRegionSecond

# TBD: compute statistics for the input and difference images created here.

# Write the extracted subimage and its difference image to disk
# For now, use 'gdalinfo' to check image statistics

writeRaster(rEdgeRegionFirst,filename="EdgeRegionFirst.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
writeRaster(rEdgeRegionAsterCgiar,filename="EdgeRegionSecond.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
writeRaster(rDelta,filename="EdgeRegionDelta.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)

[formatting edit by regetz, 06-May-2011]

I am (finally) working through the 'image difference' analysis that Jim sketched out earlier this week.
And the first thing that I did was, as he recommends, resample the CDEM image into the 'extent space' of
the ASTER/CGIAR mosaic. As expected, the statistics for the 'before' and 'after' images are very close
(note that as expected) the image sizes and resolutions did not change during the resampling):

gdalinfo -stats output - CDEM before:

reeves@eos:/data/project/organisms/rcr/ValidateBoundary$ gdalinfo -stats CanadaDemMosTuesdayClip.tif | more
Driver: GTiff/GeoTIFF
Files: CanadaDemMosTuesdayClip.tif
       CanadaDemMosTuesdayClip.aux
       CanadaDemMosTuesdayClip.rrd
Size is 41582, 3738
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-135.987917500000009,61.579680199999999)
Pixel Size = (0.000833300000000,-0.000833300000000)
Metadata:
  TIFFTAG_SOFTWARE=IMAGINE TIFF Support
Copyright 1991 - 1999 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: etif.c $ $Revision: 1.10.1.9.1.9.2.11 $ $Date: 2004/09/15 18:42:01EDT $
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-135.9879175,  61.5796802) (135d59'16.50"W, 61d34'46.85"N)
Lower Left  (-135.9879175,  58.4648048) (135d59'16.50"W, 58d27'53.30"N)
Upper Right (-101.3376369,  61.5796802) (101d20'15.49"W, 61d34'46.85"N)
Lower Right (-101.3376369,  58.4648048) (101d20'15.49"W, 58d27'53.30"N)
Center      (-118.6627772,  60.0222425) (118d39'46.00"W, 60d 1'20.07"N)
Band 1 Block=64x64 Type=Int16, ColorInterp=Gray
  Min=0.000 Max=2799.000 
  Minimum=0.000, Maximum=2799.000, Mean=613.351, StdDev=420.457
  NoData Value=-9999
  Overviews: 10396x935, 5198x468, 2599x234, 1300x117, 650x59
  Metadata:
    STATISTICS_MINIMUM=0
    STATISTICS_MAXIMUM=2799
    STATISTICS_MEAN=613.35056498368
    STATISTICS_MEDIAN=440
    STATISTICS_MODE=0
    STATISTICS_STDDEV=420.45694241697
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_EXCLUDEDVALUES=-9999
    LAYER_TYPE=athematic

gdalinfo -stats output - Resampled CDEM:

reeves@eos:/data/project/organisms/rcr/ValidateBoundary$ gdalinfo -stats CDemMosTuesdayClipMergeSpace.tif
Driver: GTiff/GeoTIFF
Files: CDemMosTuesdayClipMergeSpace.tif
       CDemMosTuesdayClipMergeSpace.aux
       CDemMosTuesdayClipMergeSpace.rrd
       CDemMosTuesdayClipMergeSpace.tif.aux.xml
Size is 41582, 3738
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-135.988060499999989,61.579322400000002)
Pixel Size = (0.000833300000000,-0.000833300000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-135.9880605,  61.5793224) (135d59'17.02"W, 61d34'45.56"N)
Lower Left  (-135.9880605,  58.4644470) (135d59'17.02"W, 58d27'52.01"N)
Upper Right (-101.3377799,  61.5793224) (101d20'16.01"W, 61d34'45.56"N)
Lower Right (-101.3377799,  58.4644470) (101d20'16.01"W, 58d27'52.01"N)
Center      (-118.6629202,  60.0218847) (118d39'46.51"W, 60d 1'18.78"N)
Band 1 Block=41582x1 Type=Int16, ColorInterp=Gray
  Min=0.000 Max=2796.000 
  Minimum=0.000, Maximum=2796.000, Mean=613.354, StdDev=420.415
  Overviews: 10396x935, 5198x468, 2599x234, 1300x117, 650x59
  Metadata:
    LAYER_TYPE=athematic
    STATISTICS_MINIMUM=0
    STATISTICS_MAXIMUM=2796
    STATISTICS_MEAN=613.35357465632
    STATISTICS_STDDEV=420.41456798272

For reference, extent of merged ASTER/SRTM image:

Corner Coordinates:
Upper Left  (-135.9880605,  61.5793224) (135d59'17.02"W, 61d34'45.56"N)
Lower Left  (-135.9880605,  58.4644470) (135d59'17.02"W, 58d27'52.01"N)
Upper Right (-101.3377799,  61.5793224) (101d20'16.01"W, 61d34'45.56"N)
Lower Right (-101.3377799,  58.4644470) (101d20'16.01"W, 58d27'52.01"N)
Center      (-118.6629202,  60.0218847) (118d39'46.51"W, 60d 1'18.78"N)

This file: CDemMosTuesdayClipMergeSpace.tif - will be the basis for the difference calculations..

[formatting edit by regetz, 06-May-2011]

Was:

gdalwarp -te -135.9880605 58.4644470 -101.3377799 61.5793224 -r bilinear CanadaDemMosTuesdayClip.tif
CDemMosTuesdayClipMergeSpace.tif

The -te coordinates are the from the gdalinfo command output for the merged Aster/CGIAR mosaic.

A note: the CDEM and merged images that I am using for this were created with a 'clip' operation,
using ArcMap GIS to create a clipping polygon covering the largest area common to both images,
and the Raster Extract By Mask command. Of course, even though clipped to the same polygon, the
origin/extent of the two images differ slightly.

#
# Usage: 
# > setwd("/data/project/organisms/rcr/ValidateBoundary/")
# > source("SampleDemDiffCols.r")
# > SampleDemDiffCols()  OR simply clip out the script between the outermost {  }

##############################################################################################
#
# SampleDemDiffCols
# 
# R script generates (or reads in) the CDEM / ASTER-SRTM mosaic 'difference image',
# and for randomly-selected one column / two row subimages, computes and saves tbe 
# difference between the pixels in the pair to create a distribution of differences.
# Three distributions created: North (ASTER CDEM), South i(SRTM/CGIAR), and 
# Border (boundary between ASTER and SRTM), Summary statistics are created for each 
# distribution.
#
# Author: Rick Reeves, NCEAS
# April 29, 2011
##############################################################################################
#
SampleDemDiffCols <- function()
{
require(raster)
require(rgdal)

#inputFirstRaster <- raster(sFirstImageName)
#inputSecondRaster <- raster(sSecondImageName)

inputFirstRaster <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")
inputSecondRaster <- raster("/data/project/organisms/rcr/ValidateBoundary/CDemMosTuesdayClipMergeSpace.tif")
#
# Difference image for entire merged image takes a while to create, 
# so we created it once, now read it back in.
#
rDeltaWhole <- raster("/data/project/organisms/rcr/ValidateBoundary/DeltaEntireImage.tif")
rDeltaWhole@data@values <-getValues(rDeltaWhole)

# Create extent objects used to extract raster subimges. 
# The object will be centered along the 60 degree North latitude line,
# and have varying depths (number of rows).
# The Western Canada study area runs from -135 (west) to -100 (west) longitude,
# and 55.0 to 64.00 degrees (north) latitude. 
# the ASTER and SRTM/CGIAR image components are merged at the 60 Deg N Latitude line.

#eTestAreaExtent <- extent(-135.2,-100.2, 59.997,60.001) # Creates a 5 row subimage
eTestAreaExtent <- extent(-135.2,-100.2, 59.995,60.005) # Creates a 12 row subimage

# Extract a sub image corresponding to the selected extent.
# Two different alternatives:
# The extract() function returns a vector of cell values, 
# the crop() function returns a complete raster* object.

vEdgeRegionFirst <- extract(inputFirstRaster,eTestAreaExtent)
rEdgeRegionFirst <- crop(inputFirstRaster,eTestAreaExtent)
vEdgeRegionSecond <- extract(inputSecondRaster,eTestAreaExtent)
rEdgeRegionSecond <- crop(inputSecondRaster,eTestAreaExtent)

# Important: In order for the image subtraction to work, the extents
#            of the two images must be IDENTICAL. I used ArcMap GIS Raster Crop By Mask
#            to create subimages with identical extents. 

# Compute the difference image  for the entire study area, and for the region along
# the boundary (narrow, maybe 10 pixels either side)

rDeltaEdge <- rEdgeRegionFirst - rEdgeRegionSecond

# Create this image one time, read it in thereafter.
#rDeltaWhole <- inputFirstRaster - inputSecondRaster
#writeRaster(rDeltaWhole,filename="DeltaEntireImage.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)

# Using the large difference image, compute subimagee statistics for areas
# North (ASTER) and South (SRTM) of the boundary. These give us an idea 
# re: differences between ASTER and CDEM and CGIAR/SRT and CDEM
# what is raster package way of using subscripts to extract? 

# Now, the interesting part: using the boundary difference image, randomly select 
# one-degree N-S strips throughout the image, and compare adjacent pixel pairs 
# above and below the boundary with pixel pairs straddling the boundary. Subtract  
# the pairs, save the collection of (absolute value) of the differences in a vector,
# so that we have a population of differences above, below, and straddling the boundary 
# line. Compare the populations.

# get a vector of random column index numbers, constrained by column dimension of image
# Loop three times, sampling pixel pairs from above, below, across the border

nColsToGet <- 1000
iDiffVecNorth <- vector(mode="integer",length=nColsToGet)
iDiffVecBorder <- vector(mode="integer",length=nColsToGet)
iDiffVecSouth <- vector(mode="integer",length=nColsToGet)

#colsToGet <-sample(1:50,nColsToGet)

# Note: initially, sample the same columns in all regions to get a profile.
#       other 'sample()' calls can be commented out to sample differenct
#       coluns in each 'region'.

# iDiffVecxxxx is a population of differences between adjacent cell pairs. 
# Compute iDiffVecNorth/Border/South on either side of border, and across it. 
# note that North and South samples taken from larger difference image for 
# entire mosaic (sub) image; iDiffBorder taken from the edge region extracted
# from the center of the lerger image.

# Remember, we are sampling a PAIR of pixels (same column from two adjacent rows)

colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet)
message("North Sample")
#browser()
# debug
#nColsToGet <- 2
#colsToGet <- c(20,100)
iFirstRow <- 300
iCtr = 1
for (iNextCol in colsToGet)
{
  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
  neighborCells <- rDeltaWhole@data@values[rColVec]
  iDiffVecNorth[iCtr] <- neighborCells[2] - neighborCells[1]
  iCtr = iCtr + 1
}
#
message("Border Sample - different columns")
#browser()
colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet)
iFirstRow <- 6 # straddle the border of 12 row center section 
iCtr = 1
for (iNextCol in colsToGet)
{
  rColVec <- cellFromRowCol(rDeltaEdge,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
  neighborCells <- rDeltaEdge@data@values[rColVec]
  iDiffVecBorder[iCtr] <- neighborCells[2] - neighborCells[1]
  iCtr = iCtr + 1 
}
#
message("South Sample - different columns")
#browser()
colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet)
iFirstRow <- 3600
iCtr = 1
for (iNextCol in colsToGet)
{
  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
  neighborCells <- rDeltaWhole@data@values[rColVec]
  iDiffVecSouth[iCtr] <- neighborCells[2] - neighborCells[1]
  iCtr = iCtr + 1 
}
# Compute iDiffVecs on either side of border, and across it. 
message("Check the cell difference vectors...")
#browser()

# summary stats for each population

sNorthSum <- sprintf("ASTER sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
                     min(iDiffVecNorth,na.rm=TRUE),median(iDiffVecNorth,na.rm=TRUE),mean(iDiffVecNorth,na.rm=TRUE),
                     max(iDiffVecNorth,na.rm=TRUE),var(iDiffVecNorth,na.rm=TRUE),sd(iDiffVecNorth,na.rm=TRUE))

sBorderSum <- sprintf("Border sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
                     min(iDiffVecBorder,na.rm=TRUE),median(iDiffVecBorder,na.rm=TRUE),mean(iDiffVecBorder,na.rm=TRUE),
                     max(iDiffVecBorder,na.rm=TRUE),var(iDiffVecBorder,na.rm=TRUE),sd(iDiffVecBorder,na.rm=TRUE))

sSouthSum <- sprintf("STRM sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
                     min(iDiffVecSouth,na.rm=TRUE),median(iDiffVecSouth,na.rm=TRUE),mean(iDiffVecSouth,na.rm=TRUE),
                     max(iDiffVecSouth,na.rm=TRUE),var(iDiffVecSouth,na.rm=TRUE),sd(iDiffVecSouth,na.rm=TRUE))
#
message(sprintf("statistics for %d N/S adjacent pixel pairs from three mosaic image regions:",nColsToGet))
message(sNorthSum)
message(sBorderSum)
message(sSouthSum)

message("hit key to write output images...")
browser()

# Write the extracted subimage and its difference image to disk
# For now, use 'gdalinfo' to check image statistics

writeRaster(rEdgeRegionFirst,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionFirstDC.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
writeRaster(rEdgeRegionSecond,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionSecondDC.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
writeRaster(rDeltaEdge,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionDeltaDC.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
}

Rick Reeves wrote:

Preliminary visual review of these files, using ArcMap GIS zoom and variable transparency tools,
shows a smooth transition between the SRTM and ASTER components.

This is not what I recall when you, Jim, and I met last Wed or so-- rather by zooming closely in to one small area, we observed what appeared to be a very distinct boundary artifact. This clearly recommended that more comprehensive and quantitative treatment is needed to determine the extent to which boundary issues were problematic, and how to develop the best strategy for integrating these two coverages.
mps

I retrieved the Wilson and Gallant text from UCSB Library, and reviewed several of the chapters.
At the very least, the first five chapters provide very good information re: the nature of the
data source(s) that we are using. Also some insights into solving the '60 degree' boundary issue.

The book is in our office, ready for your perusal.

Looking forward to discussing these with the project team.

< Rob's request is at the bottom of this posting..>
Rob:

Here is a 'distilled' set of image statistics. I will also post in Redmine, but thought to also send to you.
I will generate image histograms using R and send to you, along with R code. All of these images
are on the /Jupiter server under the subfolders 'rcr/AsterCgiarMerge' and 'rcr/ValidateBoundary'

Three groups of results.....

Rick

*********************************************************************************
1) Statistics for mosaic image components North (ASTER) and South (CGIAR / SRTM) of 60 Deg N Lat:

ASTER image component (mergeCgiarAsterBdyASTER_BL.tif)

Size is 48002, 2400
Pixel Size = (0.000833300000000,-0.000833300000000)
  Minimum=-93.000, Maximum=6366.000, Mean=665.185, StdDev=499.388
  Metadata:
    STATISTICS_MINIMUM=-93
    STATISTICS_MAXIMUM=6366
    STATISTICS_MEAN=665.18477873318
    STATISTICS_STDDEV=499.38831312001

CGIAR / SRTM image component (mergeCgiarAsterBdySRTM_BL.tif) (different (of course) and larger area than ASTER)
Size is 60002, 6000
Pixel Size = (0.000833300000000,-0.000833300000000)
 Minimum=-329.000, Maximum=4324.000, Mean=750.596, StdDev=450.075
 NoData Value=-9999
 Metadata:
  STATISTICS_MINIMUM=-329
  STATISTICS_MAXIMUM=4324
  STATISTICS_MEAN=750.59572054153
  STATISTICS_STDDEV=450.07531631289

These two files, mosaiced together (bilinear resampling) (mergeCgiarAsterBdyFinalBL.tif)
Size is 60002, 8400
Pixel Size = (0.000833300000000,-0.000833300000000)
 Minimum=-329.000, Maximum=6366.000, Mean=720.989, StdDev=469.431
 STATISTICS_MINIMUM=-329
 STATISTICS_MAXIMUM=6366
 STATISTICS_MEAN=720.98934529784
 STATISTICS_STDDEV=469.431391487922

2) Next, for a subiimage focusing on the boundary region, for which CDEM imagery is available:

Merged ASTER/SRTM Image: mergeCgiarAsterBdyTuesdayClip.tif
Size is 41582, 3738
Origin = (-135.988060499999989,61.579322388888890)
Pixel Size = (0.000833300000000,-0.000833300000000)
  Minimum=-74.000, Maximum=3851.000, Mean=623.511, StdDev=419.871
  STATISTICS_MINIMUM=-74
  STATISTICS_MAXIMUM=3851
  STATISTICS_MEAN=623.51100299321
  STATISTICS_STDDEV=419.87071089235

CDEM Image: CDemMosTuesdayClipMergeSpace.tif
Size is 41582, 3738
Origin = (-135.988060499999989,61.579322400000002)
Pixel Size = (0.000833300000000,-0.000833300000000)

Minimum=0.000, Maximum=2796.000, Mean=613.354, StdDev=420.415
  STATISTICS_MINIMUM=0
  STATISTICS_MAXIMUM=2796
  STATISTICS_MEAN=613.35357465632
  STATISTICS_STDDEV=420.41456798272

Difference Image: DeltaEntireImage.tif  (Merged ASTER/SRTM - CDEM image)
Size is 41582, 3738
Origin = (-135.988060499999989,61.579322388999998)
Pixel Size = (0.000833300000000,-0.000833300000000)
Minimum=-1603.000, Maximum=2770.000, Mean=9.047, StdDev=132.886
NoData Value=-32768
Metadata:
  STATISTICS_MINIMUM=-1603
  STATISTICS_MAXIMUM=2770
  STATISTICS_MEAN=9.0470816079223
  STATISTICS_STDDEV=132.8857558074

3) Statistics for population of randomly sampled pixel pairs (capturing North to South pixel-to-pixel change,
   as would occur along the mosaic component boundary) from the DIFFERENCE Image, computed by subtracting
   the CDEM image from the ASTER/SRTM mosaic image:

   The mean of the difference between adjacent pixels at the boundary (11.8) is close to the
   mean of the difference between the larger ASTER/SRTM image and the corresponding CDEM image (9.05).

statistics from sampling  1000 N/S adjacent pixel pairs from three mosaic image regions: 
ASTER sample summary: Min: -64.000000 / Median: 0 / Mean: 0.022000 / Max: 114.000000 / Variance: 112.043560 sDev: 10.585063
Border sample summary ***: Min: -60.000000 / Median: 11 / Mean: 11.802254 / Max: 97.000000 / Variance: 188.880856 sDev: 13.743393 
STRM sample summary: Min: -53.000000 / Median: 0 / Mean: 0.555668 / Max: 81.000000 / Variance: 96.346442 sDev: 9.815622

statistics from sampling  20000 N/S adjacent pixel pairs from three mosaic image regions:

statistics for 20000 N/S adjacent pixel pairs from three mosaic image regions:
ASTER sample summary: Min: -95.000000 / Median: 0 / Mean: -0.110150 / Max: 169.000000 / Variance: 126.035519 sDev: 11.226554
Border sample summary ***: Min: -92.000000 / Median: 11 / Mean: 11.763342 / Max: 270.000000 / Variance: 214.449517 sDev: 14.644095
STRM sample summary: Min: -83.000000 / Median: 0 / Mean: 0.397250 / Max: 115.000000 / Variance: 74.741780 sDev: 8.645333

*************************************************************************************

On 5/4/2011 12:02 PM, Robert Guralnick wrote:
> Rick --- Just send the most standard measures... not as important as means and SD --- or just show a histogram!  Visual is good.
> -r
>
> On Wed, May 4, 2011 at 12:50 PM, Rick Reeves <reeves@nceas.ucsb.edu> wrote:
>
>     Hi Rob:
>
>     I have already computed the standard summary statistics for these images, can send them now.
>
>     Could you (or Jim) be a bit more specific re: which measures of skew/kurtosos you prefer to see?
>
>     Thanks,
>     Rick
>
>
>     On 5/4/2011 11:03 AM, Robert Guralnick wrote:
>
>
>          Hey Rick and Jim --- I have gone through the Redmine conversations and generally am happy with progress.  I basically follow all the steps you have done, but I would also like a VERY SIMPLE, succinct summary of the amount of error that you are finding between ASTER and SRTM and the fused ASTER/SRTM and CDEM --- I found it hard to pull out those numbers from the outputs provided ... maybe just means, SDs, skew, kurtosis?  Or some way to get a handle on the distributions for those differentials between rasters?   I might appreciate it, and I am guessing Brian and Walter might as well.  Does this make sense to you?
>         Best, Rob
>
>
>

I have created and checked in a readme file which documents the steps I used to
generate one of the first image mosaics that I am using for exploring the DEM
data set components.
Here is the file: ReadmeCreateBoundaryMosaic.txt, located along with the scripts
that it documents at: /data/project/organisms/rcr/AsterCgiarMerge/.

Note that the specific parameters passed to 'gdalwarp' will probably change
as I adapt the process to the lessons learned working with these data sets
and creating the first image mosaics. But these scripts and this ReadMe file
demonstrate the image mosaic creation process.

May 5, 2011 / Rick Reeves
ReadMe file: Shell scripts to create the mosaic (3 arcsecond resolution)
of CGIAR/SRTM and ASTER GDEM imagery spanning the 60 Degree North Latitude
boundary within western Canada

Three shell scripts used:

1) "mosaicCgiarSrtmBdySRTM.sh" : creates mosaic from CGIAR/SRTM 5 degree tiles
    (3 arcsecond resolution) converted to .tif format from original ArcMap ASCII
    format.

2) "mosaicCgiarSrtmMBdyAster.sh" : creates mosaic from ASTER GDEM 1 degree tiles
   (1 arcsecond resolution) converted to .tif format from original ArcMap ASCII
   format. Note: incoming image tiles are resampled to 3 arcsecond resolution
   to create outgoing image mosaic,

3) "mosaicSrtmAsterPartsBdy.sh" : creates mosaic from the outputs from scripts
    1) and 2) (3 arcsecond resolution):
                                            mergeCgiarAsterBdySRTM_BL.tif
                                            mergeCgiarAsterBdyASTER_BL.tif

To create the ASTER / CGIAR boundary mosaic:

1) make current directory '/data/project/organisms/rcr/AsterCgiarMerge/'
2) run the script "mosaicCgiarSrtmBdySRTM.sh" 
3) run the script "mosaicSrtmMBdyAster.sh" 
4) run the script "mosaicSrtmAsterPartsBdy.sh" 

Note: These initial scripts do NOT specify a new image extent for any of the mosaic components.
      (gdalwarp parameter '-te xmin ymin xmax ymax' used to specify the extent)
      So for these initial *test* runs, the output image extent is derived from the input files.

      In 'production' runs, the extent of output mosaic tile compoents will be specified using -te;
      we may also set other 'gdalwarp' parameters to properly establish other output image parameters.

Two points regarding DEM raster extents and alignments:

• At least for SRTM, the pixel centers are lined up along degree boundaries. This means that the pixel borders (and thus pixel corner coordinates) are offset from the exact degree boundaries by 1/2 cell resolution. This is not a problem per se, but I mention it as a partial explanation of something to be aware of. I didn't check whether this is the case for ASTER because I wasn't sure exactly which ASTER files you are using as inputs.

• When you create an output raster by specifying the x and y cell resolution, I imagine it could cause problems because of rounding error: there is no way to specify 3/3600 (for the 90m) or 30/3600 (for the 1km) exactly in decimal form. Thus, I think it will make more sense to specify the number of pixels in the x and y direction. Specifically, this would mean using the -ts flag rather than the -tr flag in the gdal commands. Of course, you'll first have to calculate the correct number of pixels for a given output extent, but this should be easy enough.

Jim, I agree. You can in fact see the 1/2 pixel 'offset' in the screenshots that I captured
of the magnified ArcMap GIS screen.

So I will try the -ts flag on the boundary validation test case.

Workflow description: Performing the Aster/SRTM mosiac 'boundary analysis'
Rick Reeves / May 9, 2011

Objective of this analysis: explore and compare the statistical attributes of the
two Digital Elevation Model products - ASTER GDEM (north) and CGIAR SRTM (south) 
with the designated 'baseline' DEM product: Canada DEM (CDEM) (spans border) -
In proximity to the 60 degree North Latitude line, where ASTER and SRTM data sets
are merged to form the Fused Digital Terrain Layer product.

Steps, including scripts used at each stage:

1) Create three image mosaics adjacent to and spanning 60 degree North Latitude region in Western Canada: 
      CGIAR / SRTM         (north of boundary)
      ASTER GDEM           (south of boundary)
      CDEM                 (spanning boundary)
      CGIAR / ASTER mosaic (spanning boundary)

      Shell scripts used to create these image mosaics
      (in folder "/data/project/organisms/rcr/AsterCgiarMerge") 
         scripts typically executed in the following order -

         1) mosaicCgiarSrtmBdySRTM_Fix.sh  - creates SRTM image mosaic
         2) mosaicCgiarSrtmBdyAster_Fix.sh - creates Aster image mosaic
         3) mosaicSrtmAsterPartsBdyFix.sh  - creates merged SRTM/Aster image mosaic

      Creates baseline CDEM image for computing 'difference images'
         4) mosaicCgiarAsterBdyCDEM_Fix.sh

       Note: Input and output .tif image filenames are shown in the respective .sh files...

2) Displayed all of the output files produced by these scripts using ArcMapGIS interactive viewing tools
     - Confirm that the boundary area subimages are correctly located, have correct # rows, etc.

3) Developed R scripts to explore the DEM data sources adjacent to 60 degree North latitude boundary
   (in folder "/data/project/organisms/rcr/ValidateBoundary")
     Scripts typically executed in the following order...

     extractDemSubimages.r     - Extract subimages from DEM data sources for Boundary Analysis
     SampleDemDiffCols.r       - Extracts and generates statistics for randomly selected image pixel pairs
     makeImagePairPlotTable.r  - Randomly samples image columns to generate data for histograms and plots

     Note: Input and output image file nasmes are shown in the respective R scripts.

4) Use R and GDAL to generate plots and image statistics for 'images of interest' along boundary area. 
   Post results to Redmine, Iplant Wiki, etc.

5/10/2011
Rick R -

To (finally) completely respond to Rob G's request:
(Work was delayed because of unexpected problems with R Raster sub image extraction ('different origin' error)
I had to use gdalwarp to extract the subimages, then read them in to R to create 'delta' images

Added R script to above workflow: 'OnlyImageDifferences.r'

Image statistics for Aster, SRTM, CDEM, and 'Difference Images' (Aster-CDEM and SRTM-CDEM) 
for 60 Deg North Latitudeboundary region

First, for comparison, the BIG images - 

**************************** mergeCgiarAsterBdyASTER_BLEven.tif                ASTER image has higher average value

Minimum=-106.000, Maximum=6384.000, Mean=665.178, StdDev=499.372
Overviews: 12001x600, 6001x300, 3001x150, 1501x75, 751x38
Metadata:
  STATISTICS_MINIMUM=-106
  STATISTICS_MAXIMUM=6384
  STATISTICS_MEAN=665.17825516819
  STATISTICS_MEDIAN=9.9582889458621e-113
  STATISTICS_MODE=1.7802139099189e-306
  STATISTICS_STDDEV=499.37215934449
  LAYER_TYPE=athematic

**************************** mergeCgiarAsterBdySRTM_BLEven.tif            SRTM image has the smallest average value

 Min=-220.000 Max=4326.000 
Minimum=-220.000, Maximum=4326.000, Mean=620.707, StdDev=467.307
NoData Value=-9999
Overviews: 12001x1500, 6001x750, 3001x375, 1501x188, 751x94, 376x47
Metadata:
  STATISTICS_MINIMUM=-220
  STATISTICS_MAXIMUM=4326
  STATISTICS_MEAN=620.70715018819
  STATISTICS_MEDIAN=4.1631777282017e-199
  STATISTICS_MODE=7.5659641201825e-307
  STATISTICS_STDDEV=467.30736407449
  LAYER_TYPE=athematic

***************************** CDEMMosCgiarAsterBdy_BLEven.tif: CDEM image      

  Minimum=0.000, Maximum=2661.000, Mean=636.938, StdDev=423.485
  NoData Value=-9999
  Overviews: 9001x900, 4501x450, 2251x225, 1126x113, 563x57
  Metadata:
  STATISTICS_MINIMUM=0
  STATISTICS_MAXIMUM=2661
  STATISTICS_MEAN=636.93846943295
  STATISTICS_STDDEV=423.48490257975

Next, the data subimages surrounding the 60 degree North Latitude Boundary

SRTM (SOUTH of boundary) Images: 

****** Files: EdgeSRTMOnly24row.tif                         Along the border, the SRTM image has a higher mean value than ASTER

Size is 42002, 24gdalinfo -sta
Coordinate System is:

  Minimum=0.000, Maximum=2039.000, Mean=587.567, StdDev=321.291
  NoData Value=-32768
  Metadata:
    STATISTICS_MINIMUM=0
    STATISTICS_MAXIMUM=2039
    STATISTICS_MEAN=587.56716743647
    STATISTICS_STDDEV=321.29063053814
    STATISTICS_MEDIAN=1.1605833708736e-198
    STATISTICS_MODE=2.2252055940782e-306
    LAYER_TYPE=athematic

*********************     Files: EdgeSrtmCDEM24row.tif      CEDM image mean higher than SRTM image mean
Size is 2002, 24

   Minimum=168.000, Maximum=2026.000, Mean=609.157, StdDev=327.189
   NoData Value=-9999
   Metadata:
  STATISTICS_MINIMUM=168
  STATISTICS_MAXIMUM=2026
  STATISTICS_MEAN=609.15666462969
  STATISTICS_STDDEV=327.18919674128

************* Files: EdgeRegionSRTMOnlyDelta24row.tif       SRTM / CDEM Difference image
Size is 40636, 24

    Minimum=-1639.000, Maximum=180.000, Mean=4.816, StdDev=21.494
    NoData Value=-32768
    Metadata:
      STATISTICS_MINIMUM=-1639
      STATISTICS_MAXIMUM=180
      STATISTICS_MEAN=4.8156474820144
      STATISTICS_STDDEV=21.494248970983

ASTER (NORTH of Border) Image 

****** Files: EdgeAsterOnly24row.tif
Size is 42002, 24

  Minimum=119.000, Maximum=2078.000, Mean=571.326, StdDev=311.359
  NoData Value=-32768
  Metadata:
    STATISTICS_MINIMUM=119
    STATISTICS_MAXIMUM=2078
    STATISTICS_MEAN=571.32635053093
    STATISTICS_STDDEV=311.35914650438
    STATISTICS_MEDIAN=1.1610171431127e-198
    STATISTICS_MODE=9.3460979006185e-307
    LAYER_TYPE=athematic

************************  Files: EdgeAsterCDEM24row.tif  CDEM image mean higher than ASTER image mean
Size is 42002, 24

 Minimum=166.000, Maximum=1988.000, Mean=605.707, StdDev=323.313
NoData Value=-9999
Metadata:
  STATISTICS_MINIMUM=166
  STATISTICS_MAXIMUM=1988
  STATISTICS_MEAN=605.70689703064
  STATISTICS_STDDEV=323.3128648326

********************* Files: EdgeRegionAsterOnlyDelta24row.tif
Size is 40636, 24

   Minimum=-307.000, Maximum=235.000, Mean=-9.008, StdDev=29.569
   NoData Value=-32768
   Metadata:
     STATISTICS_MINIMUM=-307
     STATISTICS_MAXIMUM=235
     STATISTICS_MEAN=-9.0081548660685
     STATISTICS_STDDEV=29.568978110945

To compare performance of the candidate solution techniques, I will compute three derived terrain metrics -
Slope, Aspect, and Flow Direction - using the baseline CDEM layer. Then, I will compute them using the other seven mosaic components. Then, I will generate difference images, truth tables, RMSE, and Cross-correlation metrics, using the baseline derived components as 'truth', and the mosaic components as 'estimates' or 'simulated' values.

I will use R and MATLAB (if necessary) to generate the metrics.