[R-sig-Geo] Modis tiles with different resolutions
Matteo Mattiuzzi
matteo.mattiuzzi at boku.ac.at
Mon May 13 12:42:57 CEST 2013
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
More information about the R-sig-Geo
mailing list