projections and grid resolution

[Steve Cumming]

In the LandWEB application, which of the input data sources is taken to provide the projection for the spatial data?

If it is the ecoregion shapefile, how is the resolution of the raster maps preserved? or does it not matter?

In many of my models, we need to determine cell area, and I have been doing this from the raster resolution of a raster derived from the spTransformed vegMap (as in scfmCrop)

This worked when the studyArea shapefiles were taken from the vegMap itself (eg by randomPolygon), not when the studyArea shapefile is from some external source, in geographic coordinates. 

I am not sure what is to be done other than to abandon the idea that the input shapefile should determine the projections.

Has anyone any thoughts?

sc

[Steve Cumming]

So i have a general work-around which is to do this:

cellSize <- mean(raster::values(area(sim$flammableMap,na.rm=TRUE)),

                   na.rm=TRUE)*100 #area returns km^2!

a little slow, but has the merit of being more precise. 

The thing is, I am pretty sure it did not used to be necessary. Something has changed, I suppose in my *(&#&@?&# code, but it tracks the same as the original version written by EBJM a couple years ago, so I don't know wtf.

sc

[Alex Chubaty]

Without some sample code, or better yet, a reproducible example, I can't help you out. Generally, reprojecting an existing raster to match the projection of a shapefile keeps the raster resolution. See the ?projectRaster help file for details.

You may also want to use gdalUtils::gdalwarp for more control directly using GDAL. We also have been playing around with velox and fasterize for quick rasterization. See SpaDES.tools::fastRasterize.

[Eliot McIntire]

If your object is a raster, then:

res(rasterObj) 

gives you the resolution (vector, length 2, because you can have rasters in which the 2 dimensions are not of equal size)

In the LandWeb project, we use an input study area (called shpStudyRegionFull), which gives us the "Large Relevant Study Area", and a sub region of that (shpStudySubRegion) which we use during development, or for any time we only need a subset of the full area. Both of those are shapefiles, so as you say, do not have a resolution.  We take initial *extent* from these two objects. The extents are actually not quite exact for a raster, and these actually get updated (see next paragraph).

There are MANY rasters used throughout the LandWeb project. The first one that is loaded, called biomassMap, is taken as the "canonical" raster. We crop and mask this file with the shpStudySubRegion as the "cookie cutter". From that point on, we use biomassMap resolution ( res(sim$biomassMap) ), coordinate reference system (crs(sim$biomassMap)), and we make a slight change to the study extent and bottom-left corner to match this biomassMap raster. Because a raster has a resolution, the extent that came from the shapefiles above will not match to the 250m x 250m resolution of the raster (it will be within 250m of that) because a shapefile can have real number coordinates, and thus real number extent.

For biomassMap and all the remaining rasters we load, we use prepInputs which can specify both 

Cache(prepInputs, ...

  studyArea = sim$shpStudySubRegion # used to mask, i.e,. put NA in all pixels outside of the sim$shpStudySubRegion cutline

  rasterToMatch = sim$biomassMap # used for resolution, crs, extent

)

which will accomplish the above.

[Steve Cumming]

Yes, Alex, that was my understanding, and that is how believe the code below used to work:

 vegProjection <- crs(sim$vegMapInit) #LCC05

  

  if (is.na(crs(sim$studyArea)))            #in case it was sampled from the vegmap.

    crs(sim$studyArea) <- vegProjection

  

  simProjection<-crs(sim$studyArea)         #Ecoregion shapefile

  

  #if(ncell(sim$vegMap)>5e5) beginCluster(min(parallel::detectCores(),6))

  

  #Project the study area into each input raster, then crop and mask; 

  #Then project result back into the sim projection.

 

  studyAreaTmp <- spTransform(sim$studyArea, CRSobj =vegProjection)

  sim$vegMap <-  crop(sim$vegMapInit, studyAreaTmp)

  crs(sim$vegMap) <- vegProjection

  sim$vegMap <- mask(sim$vegMap,studyAreaTmp) #

  sim$vegMap <- projectRaster(sim$vegMap,crs=simProjection,method="ngb")

  sim$Mask <- vegMap

  sim$Mask[] <- ifelse(is.na(sim$vegMap[]), NA, 1)

  

  tmp <- getColors(sim$vegMapInit)[[1]]        # mask removes colors!

  setColors(sim$vegMap, n=length(tmp)) <- tmp  # so put them back.

  

[Steve Cumming]

I know how to find the resolution, what I cant find is what is changing it.

How reliable is prepInputs? Because I see the advantage of using a shapefile for one purpose (crop and mask) and a master-raster to specifiy resolution etc,

sc

[Alex Chubaty]

projectRaster has a 'to' argument, which is how you specify things like res, ext, and crs. See ?projectRaster.

You aren't passing 'to', so it doesn't know what to do for res. The solution is to pass a template raster to the 'to' argument. See example below:

library(raster)

f1 <- "~/Documents/shared/age/age.tif"

r1 <- raster(f1)

crs(r1) ## +proj=aea ...

res(r1) ## 1000 1000

f2 <- "~/Documents/shared/LandCoverOfCanada2005_V1_4/LCC2005_V1_4a.tif"

r2 <- raster(f2)

crs(r2) ## +proj=lcc

res(r2) ## 250 250

## this is how you've been doing it -- no "to" argument specified

r3 <- projectRaster(from = r2, crs = crs(r1))

crs(r3) ## +proj=aea

res(r3) ## 265 250 <-- uh-oh!

## the solution -- pass a template raster as the "to" argument

r4 <- projectRaster(from = r2, to = r1)

crs(r4) ## +proj=aea

res(r4)  ## 1000 1000

[Eliot McIntire]

As far as I can tell, we all have the same issues, give or take. That is why we made prepInputs. It is likely *more* reliable than anything you will come up with (prepInputs is the culmination of dozens of attempts at doing this exact thing).

It likely will not do everything you need, but it will get you close, then you modify the output of prepInputs.

The *MUCH* more sneaky thing about the raster package is that in the function "projectRaster", there is one argument "to" and several others including "crs". There is no combination of "the other argumenets" that will replicate what "to" does.  This is because the arguments are: res, crs, method="bilinear",  alignOnly=FALSE, over=FALSE, filename="" , and none of this is resolution. So, if you use projectRaster without using "to", then you will likely change the resolution of the output map, EVEN if you think you aren't changing it.

Solution: use prepInputs (which uses rasterToMatch internally) or use the "to" argument in projectRaster. Note: these two options use the same "to" argument internally.

[Steve Cumming]

that's what I get for slavishly copying your old code :-(

Thanks.

I asked about prepInputs because the man page is covered with skulls and crossbones.