[R-sig-Geo] Problem with NDVI difference with subset image

[Sarah Goslee]

Hi,

On my system, it's a writeRASTER() options issue. I don't use raster
very often, but I played with the options enough to come up with a
working solution.



## this doesn't work in the right way - something about the format
## used by writeRaster isn't compatible with lssub()
#Create geotiff
writeRaster(stackIm1, filename="stackIm1.tif", format="GTiff",overwrite=TRUE)
## also, lssub() creates a file but does not return anything, so no
## need to assign it to an object
stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif", centerx
= 0, centery = 0, widthx = 50, widthy = 50)

## this works, confirming that it's a writeRaster issue
stackIm1a <- readGDAL("stackIm1.tif")
writeGDAL(stackIm1a, "stackIm1a.tif")
lssub("stackIm1a.tif", "stackIm1a.sample.tif", centerx = 0, centery =
0, widthx = 50, widthy = 50)

## This also works, but note the name of the output file
writeRaster(stackIm1, filename="stackIm1b.tif",
format="GTiff",overwrite=TRUE, bylayer=TRUE)
lssub("stackIm1b_1.tif", "stackIm1b.sample.tif", centerx = 0, centery
= 0, widthx = 50, widthy = 50)

Sarah

On Sun, Jun 23, 2013 at 4:08 PM, ASANTOS <alexandresantosbr at yahoo.com.br> wrote:
> Hi Dra. Goslee,
>
>             I work in OS linux and I like very much your lssub function,
> despites your notice. I think that the problems are in linux system,
> because:
>
>> landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))
> Error in .local(.Object, ...) :
>   `/home/asantos/Documentos/Sensoriamento remoto percevejo
> bronzeador/Contornos/tackIm2.sample.tif' does not exist in the file system,
> and is not recognised as a supported dataset name.
>
> Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",
> :
>   Cannot create a RasterLayer object from this file. (file does not exist)
>
>
>              When I try to look the geotiff created (stackIm1.sample.tif)
> there are something wrong, because the pictures doesn't open (complete
> example below).
>
> Thanks for your attention,
>
> Alexandre
>
>
> require(raster)
> require(sp)
> require(rgdal)
> require(landsat)
> #
> #Create raster
> r <- raster(nc=1000, nr=1000)
> set.seed(20130622)
> stackIm1 <- stack(lapply(1, function(x) setValues(r, round(runif(ncell(r))*
> 255))))## Simulation red band
> stackIm2 <- stack(lapply(1, function(x) setValues(r, round(runif(ncell(r))*
> 255))))## Simulation nir
>
> # define projection system
> r.geo <- CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
> # geographical datum WGS84
> proj4string(stackIm1) <- r.geo
> proj4string(stackIm2) <- r.geo
> #
>
> #Create geotiff
> writeRaster(stackIm1, filename="stackIm1.tif",
> format="GTiff",overwrite=TRUE)
> writeRaster(stackIm2, filename="stackIm2.tif",
> format="GTiff",overwrite=TRUE)
> #
>
> #Subset a geotiff image 50 x 50 pixels
> stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif", centerx = 0,
> centery = 0, widthx = 50, widthy = 50)
> stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif", centerx = 0,
> centery = 0, widthx = 50, widthy = 50)
> #
> landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))
>
>
>
>
>
> Em 23/06/2013 12:59, Sarah Goslee escreveu:
>
>> Hi,
>>
>> What happens?
>>
>> Do you get an error? Where?
>>
>> What is your sessionInfo()?
>>
>> As it says in the lssub() help, this function was written for a
>> particular purpose, is only known to work on linux, and may not be
>> widely applicable.
>>
>> It's a "use at your own risk" kind of function, and you may well be
>> better off using one of the many other methods available for
>> subsetting raster data, which would save you the whole "export as
>> GeoTIFF, subset, reimport" sequence.
>>
>> Sarah
>>
>> On Sun, Jun 23, 2013 at 11:36 AM, Alexandre Santos
>> <alexandresantosbr at yahoo.com.br> wrote:
>>>
>>> Dear Members,
>>>
>>>
>>> I'm having trouble calculating NDVI difference, first the images Geotiff
>>> can not view using Ubuntu 4.12 64-bit and therefore can not get to the NDVI,
>>> follow an example:
>>>
>>> require(raster)
>>> require(sp)
>>> require(rgdal)
>>> require(landsat)
>>> #
>>> #Create raster
>>> r <- raster(nc=1000, nr=1000)
>>> set.seed(20130622)
>>> stackIm1 <- stack(lapply(1, function(x) setValues(r,
>>> round(runif(ncell(r))* 255))))## Simulation red band
>>> stackIm2 <- stack(lapply(1, function(x) setValues(r,
>>> round(runif(ncell(r))* 255))))## Simulation nir
>>>
>>> # define projection system
>>> r.geo <- CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
>>> # geographical datum WGS84
>>> proj4string(stackIm1) <- r.geo
>>> proj4string(stackIm2) <- r.geo
>>> #
>>>
>>> #Create geotiff
>>> writeRaster(stackIm1, filename="stackIm1.tif",
>>> format="GTiff",overwrite=TRUE)
>>> writeRaster(stackIm2, filename="stackIm2.tif",
>>> format="GTiff",overwrite=TRUE)
>>> #
>>>
>>> #Subset a geotiff image 50 x 50 pixels
>>> stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif", centerx =
>>> 0, centery = 0, widthx = 50, widthy = 50)
>>> stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif", centerx =
>>> 0, centery = 0, widthx = 50, widthy = 50)
>>> #
>>>
>>> #Calculate NDVI difference
>>>
>>> multi.espc<-stack(c("stackIm1.sample.tif","stackIm2.sample.tif"))
>>>
>>> band3<-raster(multi.espc,1)
>>> band4<-raster(multi.espc,2)
>>>
>>> ndvi<-(band4-band3)/(band4+band3)
>>>
>>> #
>>>

-- 
Sarah Goslee
http://www.functionaldiversity.org