How to write functions with mutiple large raster data file (MOD13Q1, 2000-2012)

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130507/8914d520/attachment.pl>

Francisco,

I think what you want to do can be accomplished with 'calc' (as was
pointed out in another thread).

library(raster)
in <- stack(list.files(pattern=*..tif$))
out <- calc(in, fun=function(x) lowess(x,f=7/n,iter=10)$y)


But to answer your question, I think you could do:

in <- stack(list.files(pattern=*..tif$))
# and, here is the trick
out <- brick(in, values=FALSE)

Also, in your loop you could do:

getValues(in, ...
writeValues(out, ....


Best, Robert

Hi everybody,

I'm working with MODIS data MOD13Q1 and the NDVI index between 2000- today.
I have been trying to do a function, but I have a problem when stack the
310 files and then I tried to use the writeStart function. I had an error
because the startFuction work with rasterLayer or RasterBrick. But, when I
tried to make the brick with the rasterStack (>brick(rasterStack)) the
operation took too much time (11 hours)

My question, is what is the most efficiently way to make a rasterBrick from
multiple raster (geoTIFF) files?

save the rasterStack and then open it with brick?

My code is something like:

out<-stack(list.files(pattern=*..tif$))
out<-brick(mod) #too much time (11 hours appox)

  >out<-writeStart(out,filename,overwrite=TRUE)
  >bs<-blockSize(out)
  >pb<-pbCreate(bs$n)
  >for (i in 1:bs$n){

    v<-getValues(out,row=bs$row[i],nrows=bs$nrows[i])
    s<-t(sapply(apply(v,1,lowess,f=7/n,iter=10),"[[","y"))
    out<-writeValues(out,s,bs$row[i])
   pbStep(pb,i)
  }

Kind Regards

Francisco Zambrano Bigiarini

from Chile, latin america

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130508/880271e7/attachment.pl>

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130508/f243e390/attachment.pl>

Dear all,

I am processing large sets of MODIS data that come to me in sinusoidal 
projection and .lan (ERDAS 7.5) format. I have to mosaic different tiles 
and here comes the problem I meet: resolutions are not the same.

Basically I'm looping through these steps:

library(raster)
library(sp)
wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
sinu <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 
+units=m +no_defs"
## MA (h22v10) tile extent in longlat
e <- extent(c(40.6187, 53.2072, -20, -10))
## transform to sinusoidal
e <- as.data.frame(t(bbox(e)))
coordinates(e) <- ~ s1 + s2
projection(e) <- wgs84
e <- spTransform(e, CRS(sinu))
e <- coordinates(e)
e <- extent(c(e[,1], e[,2]))
## First raster
r <- raster("c:/Temp/RtmpyQ3t4R/MA050107.lan")
projection(r) <- sinu
extent(r) <- e

at the end of the process:

 > ## tile MA (h22v10)
 > MAdlst[[1]]
class       : RasterLayer
dimensions  : 1200, 1200, 1440000  (nrow, ncol, ncell)
resolution  : 1318.567, 926.6254  (x, y)
extent      : 4244214, 5826494, -2223901, -1111951  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 
+b=6371007.181 +units=m +no_defs
data source : c:\Temp\RtmpyQ3t4R\MA050107.lan
names       : MA050107
values      : -32768, 32767  (min, max)

 > ## tile MB (h22v11)
 > MBdlst[[1]]
class       : RasterLayer
dimensions  : 1200, 1200, 1440000  (nrow, ncol, ncell)
resolution  : 1611.032, 926.6254  (x, y)
extent      : 4099193, 6032431, -3335852, -2223901  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 
+b=6371007.181 +units=m +no_defs
data source : c:\Temp\RtmpyQ3t4R\MB050107.lan
names       : MB050107
values      : -32768, 32767  (min, max)


Am I doing something wrong, and how can I mosaic these rasters with 
different resolutions?

All the best,

Renaud


 > sessionInfo()
R version 3.0.0 (2013-04-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252 
LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods 
   base

other attached packages:
[1] rgdal_0.8-9     maptools_0.8-23 foreign_0.8-53  lattice_0.20-15 
mapdata_2.2-2   maps_2.3-2      raster_2.1-25
[8] sp_1.0-9

loaded via a namespace (and not attached):
[1] compiler_3.0.0 fortunes_1.5-1 tools_3.0.0








Le 08/05/2013 21:45, Francisco Zambrano a ?crit :

Robert,

Sorry, now I realize that calc works in chunk if is needed. And I imagine
that calc is more efficient than do it the chunks by myself.

Thanks,

Francisco


2013/5/8 Robert J. Hijmans<r.hijmans at gmail.com>

Francisco,

I think what you want to do can be accomplished with 'calc' (as was
pointed out in another thread).

library(raster)
in<- stack(list.files(pattern=*..tif$))
out<- calc(in, fun=function(x) lowess(x,f=7/n,iter=10)$y)


But to answer your question, I think you could do:

in<- stack(list.files(pattern=*..tif$))
# and, here is the trick
out<- brick(in, values=FALSE)

Also, in your loop you could do:

getValues(in, ...
writeValues(out, ....


Best, Robert

On Tue, May 7, 2013 at 7:35 AM, Francisco Zambrano<frzambra at gmail.com>
wrote:

Hi everybody,

I'm working with MODIS data MOD13Q1 and the NDVI index between 2000-

today.

I have been trying to do a function, but I have a problem when stack the
310 files and then I tried to use the writeStart function. I had an error
because the startFuction work with rasterLayer or RasterBrick. But, when

I

tried to make the brick with the rasterStack (>brick(rasterStack)) the
operation took too much time (11 hours)

My question, is what is the most efficiently way to make a rasterBrick

from

multiple raster (geoTIFF) files?

save the rasterStack and then open it with brick?

My code is something like:

out<-stack(list.files(pattern=*..tif$))
out<-brick(mod) #too much time (11 hours appox)

   >out<-writeStart(out,filename,overwrite=TRUE)
   >bs<-blockSize(out)
   >pb<-pbCreate(bs$n)
   >for (i in 1:bs$n){

     v<-getValues(out,row=bs$row[i],nrows=bs$nrows[i])
     s<-t(sapply(apply(v,1,lowess,f=7/n,iter=10),"[[","y"))
     out<-writeValues(out,s,bs$row[i])
    pbStep(pb,i)
   }

Kind Regards

Francisco Zambrano Bigiarini

from Chile, latin america

         [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Renaud Lancelot
Directeur adjoint / Deputy director
CIRAD, UMR15, Campus International de Baillarguet TA A-DIR / B
F34398 Montpellier

EDENext Project, coordinator: http://www.edenext.eu/

Tel.  +33 4 67 59 37 17  -  Fax  +33 4 67 59 37 98
Secr. +33 4 67 59 37 37  - Cell. +33 6 77 52 08 69

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130513/525a844c/attachment.pl>

Lancelot,
The way you have formulated your question and the code you have provided
does not give enough information to help in an easy and appropriate way.
Please make sure that the code you provide is enough to reproduce the
example or, if not possible, that any information is included.


I think you have mixed up something with your extent or better with its
re-projection! 
An easier way to re-project the extent of a raster* object is:
#############
library(raster)
e <- raster(xmn=40.6187, xmx=53.2072, ymn=-20, ymx=-10)
e <- projectExtent(e,"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs")
e # this 'e' is right and not the same as yours!


# Have you tested to load your results in a GIS? You will see that it
isn't just a resolution problem but the two images are not on the right
place! 
# Find bellow the extent of the two tiles using the MODIS package in an
example run (in your case probably not usable as you are working with
already available '.lan' files).
############

library(MODIS)

# TILE h22v11

runGdal(product="MOD11A2", begin="2010001", end="2010007", tileH=22,

tileV=11, SDSstring=1, job="LST")
pixelSize        =  asIn 
resamplingType   =  near 
outProj          = +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs 
Output Directory = ~/MODIS_PROCESSED/LST

raster('~/MODIS_PROCESSED/LST/MOD11A2.A2010001.LST_Day_1km.tif' )

class       : RasterLayer 
dimensions  : 1200, 1200, 1440000  (nrow, ncol, ncell)
resolution  : 926.6254, 926.6254  (x, y)
extent      : 4447802, 5559753, -3335852, -2223901  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs 
data source : ~/MODIS_PROCESSED/LST/MOD11A2.A2010001.LST_Day_1km.tif 
names       : MOD11A2.A2010001.LST_Day_1km 
values      : 0, 65535  (min, max)


# TILE h22v10

runGdal(product="MOD11A2", begin="2010001", end="2010007", tileH=22,

tileV=10, SDSstring=1, job="LST2")
pixelSize        =  asIn 
resamplingType   =  near 
outProj          = +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs 
Output Directory = ~/MODIS_PROCESSED/LST2

raster('~/MODIS_PROCESSED/LST2/MOD11A2.A2010001.LST_Day_1km.tif' )

class       : RasterLayer 
dimensions  : 1200, 1200, 1440000  (nrow, ncol, ncell)
resolution  : 926.6254, 926.6254  (x, y)
extent      : 4447802, 5559753, -2223901, -1111951  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs 
data source : ~/MODIS_PROCESSED/LST2/MOD11A2.A2010001.LST_Day_1km.tif 
names       : MOD11A2.A2010001.LST_Day_1km 
values      : 0, 65535  (min, max)


# IHTH
# Matteo



# PS:
# both TILES + mosaicked:

runGdal(product="MOD11A2", begin="2010001", end="2010007", tileH=22,

tileV=10:11, SDSstring=1, job="LST3")


# entire Madagascar:

runGdal(product="MOD11A2", begin="2010001", end="2010007",

extent="Madagascar", SDSstring=1, job="LST4")

lancelot  05/09/13 2:07 AM >>>

Dear all,

I am processing large sets of MODIS data that come to me in sinusoidal 
projection and .lan (ERDAS 7.5) format. I have to mosaic different tiles

and here comes the problem I meet: resolutions are not the same.

Basically I'm looping through these steps:

library(raster)
library(sp)
wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
sinu <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181

+units=m +no_defs"
## MA (h22v10) tile extent in longlat
e <- extent(c(40.6187, 53.2072, -20, -10))
## transform to sinusoidal
e <- as.data.frame(t(bbox(e)))
coordinates(e) <- ~ s1 + s2
projection(e) <- wgs84
e <- spTransform(e, CRS(sinu))
e <- coordinates(e)
e <- extent(c(e[,1], e[,2]))
## First raster
r <- raster("c:/Temp/RtmpyQ3t4R/MA050107.lan")
projection(r) <- sinu
extent(r) <- e

at the end of the process:

 > ## tile MA (h22v10)
 > MAdlst[[1]]
class       : RasterLayer
dimensions  : 12extent      : 4244214, 5826494, -2223901, -1111951  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 
+b=6371007.181 +units=m +no_defs
data source : c:\Temp\RtmpyQ3t4R\MA050107.lan
names       : MA050107
values      : -32768, 32767  (min, max)

 > ## tile MB (h22v11)
 > MBdlst[[1]]
class       : RasterLayer
dimensions  : 1200, 1200, 1440000  (nrow, ncol, ncell)
resolution  : 1611.032, 926.6254  (x, y)
extent      : 4099193, 6032431, -3335852, -2223901  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 
+b=6371007.181 +units=m +no_defs
data source : c:\Temp\RtmpyQ3t4R\MB050107.lan
names       : MB050107
values      : -32768, 32767  (min, max)


Am I doing something wrong, and how can I mosaic these rasters with 
different resolutions?

All the best,

Renaud


 > sessionInfo()
R version 3.0.0 (2013-04-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252 
LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods 
   base

other attached packages:
[1] rgdal_0.8-9     maptools_0.8-23 foreign_0.8-53  lattice_0.20-15 
mapdata_2.2-2   maps_2.3-2      raster_2.1-25
[8] sp_1.0-9

loaded via a namespace (and not attached):
[1] compiler_3.0.0 fortunes_1.5-1 tools_3.0.0








Le 08/05/2013 21:45, Francisco Zambrano a ?crit :

Robert,

Sorry, now I realize that calc works in chunk if is needed. And I

imagine

that calc is more efficient than do it the chunks by myself.

Thanks,

Francisco


2013/5/8 Robert J. Hijmans

Francisco,

I think what you want to do can be accomplished with 'calc' (as was
pointed out in another thread).

library(raster)
in<- stack(list.files(pattern=*..tif$))
out<- calc(in, fun=function(x) lowess(x,f=7/n,iter=10)$y)


But to answer your question, I think you could do:

in<- stack(list.files(pattern=*..tif$))
# and, here is the trick
out<- brick(in, values=FALSE)

Also, in your loop you could do:

getValues(in, ...
writeValues(out, ....


Best, Robert

On Tue, May 7, 2013 at 7:35 AM, Francisco Zambrano
wrote:

Hi everybody,

I'm working with MODIS data MOD13Q1 and the NDVI index between 2000-

today.

I have been trying to do a function, but I have a problem when stack

the

310 files and then I tried to use the writeStart function. I had an

error

because the startFuction work with rasterLayer or RasterBrick. But,

when

I

tried to make the brick with the rasterStack (>brick(rasterStack))

the

operation took too much time (11 hours)

My question, is what is the most efficiently way to make a

rasterBrick

from

multiple raster (geoTIFF) files?

save the rasterStack and then open it with brick?

My code is something like:

out<-stack(list.files(pattern=*..tif$))
out<-brick(mod) #too much time (11 hours appox)

   >out<-writeStart(out,filename,overwrite=TRUE)
   >bs<-blockSize(out)
   >pb<-pbCreate(bs$n)
   >for (i in 1:bs$n){

     v<-getValues(out,row=bs$row[i],nrows=bs$nrows[i])
     s<-t(sapply(apply(v,1,lowess,f=7/n,iter=10),"[[","y"))
     out<-writeValues(out,s,bs$row[i])
    pbStep(pb,i)
   }

Kind Regards

Francisco Zambrano Bigiarini

from Chile, latin america

         [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Renaud Lancelot
Directeur adjoint / Deputy diTel.  +33 4 67 59 37 17  -  Fax  +33 4 67 59 37 98
Secr. +33 4 67 59 37 37  - Cell. +33 6 77 52 08 69

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo