Hi,

Try raster::resample then raster::mosaic


2013/5/9 lancelot <renaud.lancelot@cirad.fr>

> 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@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@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@r-project.org
>>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>>
>>>
>>>
>>
>>
>>
>>
>> ______________________________**_________________
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<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
>
> ______________________________**_________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>

	[[alternative HTML version deleted]]

