[R-sig-Geo] is it possible to resample a raster image?
Roger Bivand
Roger.Bivand at nhh.no
Thu Apr 19 13:27:45 CEST 2007
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
More information about the R-sig-Geo
mailing list