[R-sig-Geo] is it possible to resample a raster image?

[Roger Bivand]

On Thu, 19 Apr 2007, Vladimir Eremeev wrote:

> Dear all,
> 
> I am studying new knowledge areas now.
> 
> I am curious, if it is possible to transform a projected raster image
> to another geographic projection using R packages.

There is no warping functionality as such, so as you see below, the output
data is not in a regular grid. You need to make your "+proj=longlat" (not
"latlong" - that may reverse the understanding of the coordinates)  
SpatialGrid first, extract its coordinates(), spTransform() them to the
input projection - now a SpatialPoints object. Next, interpolate from the
input grid to the new points, and finally either project back to "longlat"
- which will be a rectangular grid, or just insert the interpolated values
as the data= argument into a SpatialGridDataFrame using the "longlat"  
GridTopology created earlier.

Using GridTopology() and coordinates(), you could get round the 
expand.grid() step. For interpolation, IDW in gstat with a tight 
neighbourhood, or maybe splines in interp or field ought to be OK, or some 
other resampling trick (k=1 nearest neighbour for categorical variables?).

Hope this helps,

Roger

> 
> For example, I have a matrix (representing an image)
> It has 304 columns by 448 rows.
> Initially it is stored in a file, having 2 bytes per pixel.
> Lots of such images can be found here
> ftp://sidads.colorado.edu/pub/DATASETS/brightness-temperatures/
> polar-stereo/ssmi/north/2006/
> 
> I import one of them in R:
> 
> tb.mtx<-matrix(readBin("ftp://sidads.colorado.edu/pub/DATASETS/brightness-temperatures/polar-stereo/ssmi/north/2006/tb_f13_20060101_v2_n19h.bin",
>                        "integer",n=304*448,size=2),
>                ncol=304,byrow=TRUE)
> 
> Then I create the SpatialGridDataFrame:
> 
> egrid<-cbind(expand.grid(x=seq(-3850000.0000000,3725000.0000000,by=25000),
>                          y=seq(-5350000.0000000,5825000,by=25000)),
>              as.vector(t(tb.mtx)[,448:1]))
> 
> coordinates(egrid)<-c("x","y")
> gridded(egrid)<-TRUE
> proj4string(egrid)<-CRS("+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +ellps=clrk66")
> 
> I have worked with these data using other tools for a long time
> already, therefore, I know about their projection and dimensions.
> These are documented in the web and also in `readme` files.
> 
> Now I would like to *use R-spatial packages* to transform egrid to another projection,
> and export the result into an external flat binary file, similar to
> the one, I have imported with readBin.
> 
> I can easily do it using other tool sets, however, now I am exploring
> R.
> Is it possible?
> 
> I tried the following, however, got errors.
> > egtr<-spTransform(egrid,CRS("+proj=latlong"))
> > gridded(egtr)<-TRUE
> 
> Error in if (any(outside)) { : missing value where TRUE/FALSE needed
> In addition: Warning messages:
> 1: grid has empty column/rows in dimension 1 in: points2grid(points, tolerance) 
> 2: NAs introduced by coercion 
> 3: grid has empty column/rows in dimension 2 in: points2grid(points, tolerance) 
> 4: NAs introduced by coercion 
> 
> > fullgrid(egtr)<-TRUE
> Error in `fullgrid<-`(`*tmp*`, value = TRUE) : 
>         fullgrid<- only works on objects of class or extending SpatialPixels
> 
> 
> Thank you!
> 
> 
> ---
> Best regards,
> Vladimir                mailto:wl at eimb.ru
> 
> 
> 
> --SevinMail--
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no