[R-sig-Geo] Modis tiles with different resolutions

lancelot renaud.lancelot at cirad.fr
Thu May 9 02:07:03 CEST 2013 

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

More information about the R-sig-Geo mailing list