Skip to content

Latest commit

 

History

History
500 lines (409 loc) · 14.1 KB

pad-us_nlcd.org

File metadata and controls

500 lines (409 loc) · 14.1 KB

pad-us_nlcd.org

1 TODO tangle out the R code and run it from the Makefile

2 initialize the session

## library( raster)
library( raster, lib.loc="~/src/R/lib/")
setOptions( progress= "text")
## setOptions( maxmemory= 1e+06)
library( plyr)
library( stringr)

overwriteRasters <- FALSE

3 process Puerto Rico to work out steps

pr <- raster( "nlcd2006/pr_landcover_wimperv_10-28-08_se5.img")
NAvalue( pr) <- 0
pr <- setMinMax( pr)

prGrid <- try( raster( "prGrid.tif"), silent= TRUE)
if( inherits( prGrid, "try-error") || overwriteRasters) {             
  prGrid <- raster( pr)
  prGrid[] <- seq( 1, ncell( prGrid))
  prGrid <-
    mask( prGrid, pr,
         filename= "prGrid.tif",
         overwrite= TRUE,
         progress= "text")
}

gridProjFunc <- function( cell) {
  cellFromXY( world,
             project( xyFromCell( prGrid, cell),
                     projection( prGrid),
                     inv= TRUE))
}  

prWorld <- try( raster( "world_5min_PuertoRico.tif"), silent= TRUE)
if( inherits( prWorld, "try-error") || overwriteRasters) {             
  prWorld <-
    calc( prGrid, gridProjFunc,
         filename= "world_5min_PuertoRico.tif",
         datatype= "INT4U",
         overwrite= TRUE,
         progress= "text")
}

prGap <- raster( "pad-us/PADUS1_2_regions/PADUS1_2_PuertoRico_GAP.tif")
prGap <- setMinMax( prGap)
NAvalue( prGap) <- 255
  
prGap <- overlay( prGap, prGrid, fun= setGapZero,
                 filename= "prGap.tif", datatype= "INT1U", progress= "text", overwrite= TRUE)
NAvalue( prGap) <- 255


prStack <- stack(prWorld, pr, prGap)
layerNames( prStack) <- c( "grid", "nlcd", "gap")

ct <- crosstab( prStack, useNA= "always", long= TRUE, responseName= "n", progress="text")

4 extend PR example for batch processing

4.1 load NLCD rasters

## regions <- c( "Alaska", "Hawaii", "CUSA", "PuertoRico")
## names( regions) <- regions
## nlcdFiles <-
##   list.files( "nlcd2006",
##              patt="^(ak|hi|nlcd2006|pr).*?img$",
##              full.names= TRUE)
## names( nlcdFiles) <- regions

## regions <- regions[ c( 4, 2, 1, 3)]

## nlcdRasters <-
##   llply( regions,
##         function ( region) {
##           r <- raster( nlcdFiles[[ region]])
##           NAvalue( r) <- 0
##           ## r <- setMinMax( r)
##           layerNames( r) <- region
##           r
##         })

regionPatterns <-
  list( PuertoRico= "pr.*?img$",
        Hawaii= "hi.*?img$",
        Alaska= "ak.*?img$",
        cUSA= "nlcd2006.*?img$")

regions <- names( regionPatterns)
names( regions) <- names( regionPatterns)
        
nlcdRasters <-
  llply( regionPatterns,
        function( patt) {
          r <- raster( list.files( "nlcd2006",
                                  patt= patt,
                                  full.names= TRUE))
          NAvalue( r) <- 0
          r
        })

4.2 calculate 5’ cell ID for each 30m pixel

Write out a 5’ raster in geographic projection where the value of each cell is its grid ID. This will be reprojected into the cooridnate space of each PAD-US/NLCD stack.

world <- raster()
res( world) <- 5/60
dataType( world) <- "INT4U"
world[ ] <- seq( 1, ncell( world))
world <- writeRaster( world, "world5min.tif",
                     datatype= "INT4U",
                     overwrite= TRUE)

4.2.1 TODO How did I write the gdal_rasterize command for the grid IDs?

I must have done it by hand. This should be tangled out and called in the Makefile.

4.3 add zeroes to GAP data for unprotected land and coastal areas

gapFiles <-
  list.files( "pad-us/PADUS1_2_regions/",
             patt= "^PADUS1_2_.*?tif$",
             full.names= TRUE)
names( gapFiles) <-
  str_match( gapFiles,
            "PADUS1_2_([^_]+)_GAP\\.tif$")[, 2]

gapRasters <-
  llply( names( regionPatterns),
        function ( region) {
          r <- raster( gapFiles[[ region]])
          NAvalue( r) <- 255
          ## r <- setMinMax( r)
          layerNames( r) <- region
          r
})
names( gapRasters) <- names( regionPatterns)

setGapZero <- function( gap, grid) {
  ifelse( is.na( gap) & !is.na( grid), 0, gap)
}

gapOverlayFunc <-
  function ( gap, nlcd) {
    fn <- sprintf( "gap%s.grd", layerNames( gap))
    if( overwriteRasters | !file.exists( fn)) {
      overlay( gap, nlcd,
              fun= setGapZero,
              filename= fn,
              datatype= "INT1U",
              overwrite= TRUE)
    } else try( raster( fn), silent= TRUE)
  }

prOverlay <- gapOverlayFunc( gapRasters[[ "PuertoRico"]],
                            nlcdRasters[[ "PuertoRico"]])

## gapOverlays <-
##   mapply( gapRasters, nlcdRasters,
##          FUN= gapOverlayFunc) 

gapOverlays <-
  llply( regions,
        function( region) {
          gapOverlayFunc( gapRasters[[ region]],
                         nlcdRasters[[ region]])
        })

4.4 create stacks and tabulate

## prStack <- stack(prWorld, pr, prGap)
## layerNames( prStack) <- c( "grid", "nlcd", "gap")

## prStack <- stack( raster( "aeaGrid5minPuertoRico.img"),
##                  nlcdRasters[[ "PuertoRico"]],
##                  prOverlay)

## prLowRes <- raster( prStack)
## res( prLowRes) <- 3000

## prStackSmall <- resample( prStack, prLowRes, method= "ngb")
## layerNames( prStackSmall) <- c( "grid", "nlcd", "gap")

## prCt <- crosstab( prStackSmall, long= TRUE, responseName= "n")

## prCt <- crosstab( prStack, long= TRUE)


aeaGridFunc <-
  function( region) {
    raster( sprintf( "aeaGrid5min%s.img", region))
  }
  
aeaGrids <- llply( regions, aeaGridFunc)
                  
gapStackFunc <-
  function( region) {
    s <- stack( aeaGrids[[ region]],
               nlcdRasters[[ region]],
               gapOverlays[[ region]])
    layerNames( s) <- c( "grid", "nlcd", "gap")
    s
  }
               
gapStacks <- llply( regions, gapStackFunc)
        
writeCrosstabs <-
  function( region) {
    fn <- sprintf( "pad-us_nlcd_%s.csv", region)
    ct <- crosstab( gapStacks[[ region]])
    write.csv( ct, row.names= FALSE, file= fn)
    fn
  }

ctFiles <- llply( regions, writeCrosstabs)

4.5 write out GRASS scripts

./create_location.sh aeaGrid5minPuertoRico.img PuertoRico grass
g.rename rast=aeaGrid5minPuertoRico.img,grid_5min
r.in.gdal input=nlcd2006/pr_landcover_wimperv_10-28-08_se5.img output=nlcd2006
r.in.gdal input=pad-us/PADUS1_2_regions/PADUS1_2_PuertoRico_GAP.tif output=gap

r.mapcalc MASK="if( nlcd2006 > 0, 1, null())"
echo grid_5min,nlcd2006,gap,n > statsPuertoRico.csv
r.stats -c input=grid_5min,nlcd2006,gap fs=, >> statsPuertoRico.csv
r.mask -r
echo grid_5min,n > gridPuertoRico.csv
r.stats -c input=grid_5min fs=, >> gridPuertoRico.csv

# echo ./create_location.sh aeaGrid5min${REGION}.img ${REGION} grass

cat <<'EOF'
# path to GRASS binaries and libraries:
export GISBASE=/usr/lib/grass64
export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib

# use process ID (PID) as lock file number:
export GIS_LOCK=$$

# path to GRASS settings file
export GISRC=./.grassrc6
EOF

cat <<EOF
g.gisenv set=LOCATION_NAME=\$GIS_LOCK
eval \$(g.gisenv)
mkdir -p \$GISDBASE/\$GIS_LOCK
g.mapset -c mapset=temp
r.in.gdal in=aeaGrid5min${REGION}.img out=grid_5min location=${REGION}
# g.gisenv set=LOCATION_NAME=${REGION}
# g.gisenv set=MAPSET=PERMANENT
g.mapset mapset=PERMANENT location=${REGION}
rm -rf \$GISDBASE/\$GIS_LOCK/temp
# g.rename rast=aeaGrid5min${REGION}.img,grid_5min
r.in.gdal input=pad-us/PADUS1_2_regions/PADUS1_2_${REGION}_GAP.tif output=gap
EOF
path <-
  list.files( "nlcd2006",
             patt= regionPatterns[[ region]],
             full.names= TRUE)
sprintf( "r.in.gdal input=%s output=nlcd2006", path)

4.5.1 TODO add ‘-N’ to r.stats for NLCD/GAP tabulation to eliminate ,,*,n record created by the mask

cat <<EOF 
r.mapcalc MASK="if( nlcd2006 > 0, 1, null())"
echo grid_5min,nlcd2006,gap,n > stats${REGION}.csv
r.stats -c input=grid_5min,nlcd2006,gap fs=, >> stats${REGION}.csv
r.mask -r
echo grid_5min,n > grid${REGION}.csv
r.stats -c input=grid_5min fs=, >> grid${REGION}.csv
EOF

cat <<'EOF'
# run GRASS' cleanup routine
$GISBASE/etc/clean_temp

# remove session tmp directory:
rm -rf /tmp/grass6-$USER-$GIS_LOCK
EOF
<<grassCreate( "PuertoRico")>>
<<rInGdalNlcd( "PuertoRico")>>
<<grassMapcalc( "PuertoRico")>>
<<grassCreate( "Hawaii")>>
<<rInGdalNlcd( "Hawaii")>>
<<grassMapcalc( "Hawaii")>>
<<grassCreate( "Alaska")>>
<<rInGdalNlcd( "Alaska")>>
<<grassMapcalc( "Alaska")>>
<<grassCreate( "cUSA")>>
<<rInGdalNlcd( "cUSA")>>
<<grassMapcalc( "cUSA")>>

4.6 aggregate the results

library( reshape)
library( Hmisc)

cells <-
  read.csv( "gridPuertoRico.csv",
           col.names= c( "cell", "n"))

stats <-
  read.csv( "statsPuertoRico.csv",
           na.strings= "*",
           col.names= c( "cell", "nlcd", "gap", "n"),
           colClasses= c("numeric", "factor", "factor", "numeric"))
## won't need this when r.stats in previous GRASS step is fixed
stats <- stats[ !is.na(stats$cell),]

## stats <- stats[ !is.na(stats$grid),]
## stats <- stats[ stats$cell != "*",]

## stats <- within( stats, gap[ is.na( gap)] <- 0)


stats <-
  within( stats,
         { levels( gap) <- c( levels( gap), "0")
           gap[ is.na( gap)] <- "0"
           gap <- combine_factor( gap, c(0,1,1,1,0))
           levels( gap) <- c( "no", "yes")
         })

stats <-
  cast( data= stats,
       formula= cell ~ gap + nlcd,
       fun.aggregate= sum,
       margins= "grand_col",
       value= "n" )
colnames( stats)[ colnames( stats) == "(all)_(all)"] <- "nlcd"

merged <-
  within( merge( stats, cells, by= "cell", all.x= TRUE),
         no_11 <- no_11 + n - nlcd)

fracs <-
  cast( within( melt( merged,
                     c( "cell", "n")),
               value <- value / n),
       formula= cell ~ variable,
       subset= variable != "nlcd",
       margins= "grand_col",
       fun.aggregate= sum)

write.csv( format.df( fracs,
                     dec= 3,
                     numeric.dollar= FALSE,
                     na.blank= TRUE),
          row.names= FALSE,
          file= "fracsPuertoRico.csv",
          quote= FALSE)
library( reshape)
library( Hmisc)

writeFracs <- function( region) {
  cells <-
    read.csv( sprintf( "grid%s.csv", region),
             col.names= c( "cell", "n"))
  stats <-
    read.csv( sprintf( "stats%s.csv", region),
             na.strings= "*",
             col.names= c( "cell", "nlcd", "gap", "n"),
             colClasses= c("numeric", "factor", "factor", "numeric"))
  ## won't need this when r.stats in previous GRASS step is fixed
  stats <- stats[ !is.na(stats$cell),]
  stats <-
    within( stats,
           { levels( gap) <- c( levels( gap), "0")
             gap[ is.na( gap)] <- "0"
             gap <- combine_factor( gap, c(0,1,1,1,0))
             levels( gap) <- c( "no", "yes")
           })
  stats <-
    cast( data= stats,
         formula= cell ~ gap + nlcd,
         fun.aggregate= sum,
         margins= "grand_col",
         value= "n" )
  colnames( stats)[ colnames( stats) == "(all)_(all)"] <- "nlcd"
  merged <-
    within( merge( stats, cells, by= "cell", all.x= TRUE),
           no_11 <- no_11 + n - nlcd)
  fracs <-
    cast( within( melt( merged,
                       c( "cell", "n")),
                 value <- value / n),
         formula= cell ~ variable,
         subset= variable != "nlcd",
         margins= "grand_col",
         fun.aggregate= sum)
  fn <- sprintf( "fracs%s.csv", region)
  write.csv( format.df( fracs,
                       dec= 3,
                       numeric.dollar= FALSE,
                       na.blank= TRUE),
            row.names= FALSE,
            file= fn,
            quote= FALSE)
  fn
}

regions <- c( "PuertoRico", "Hawaii", "Alaska", "cUSA")
names( regions) <- regions

fracFiles <- llply( regions, writeFracs)

zip( "pad-us_nlcd.zip", list.files( patt= "^fracs.*?\\csv$"))