1 Jul 2010 07:02

### Re: polygonValues (raster): Very slow

```Pierre,

This looks great also, thanks. It's going to take a while, though, as
I'm working
on linux and will have to compile WKT and probably PostGis also (don't
know if the
later comes included with the binary distro of qgis).
Why don't you talk to the ubuntugis people for a regular ubuntu
package? You could use
the unstable branch if needed:
This would make things so much easier for people using ubuntu and WKT would be
tested by many more users.

I'll let you know how this works.

Agus

2010/6/30 Pierre Racine <Pierre.Racine <at> sbf.ulaval.ca>:
> Agustin,
>
> PostGIS WKT Raster was made to do exactly this kind of operation on large datasets. Give it a try!
>
> http://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01
>
> Pierre
>
>>-----Original Message-----
>>From: r-sig-geo-bounces <at> stat.math.ethz.ch [mailto:r-sig-geo-bounces <at> stat.math.ethz.ch] On
Behalf Of
```

1 Jul 2010 10:59

### make thousands SpatialPolygons()

```Hi,
I want to make approximately 10000 square polygons (300*300) based on coordinates in a data frame, and make
a .shp-file out of these.

The .pdf "S Classes and Methods for Spatial Data: the sp Package" by Pebesma and Bivand contains an
instruction on how to do this, and it works very well, e.g.

i=1
Sr1 = Polygon(cbind(c(dat[i,"X"]-150,dat[i,"X"] - 150,dat[i,"X"]+150,dat[i,"X"]+150,dat[i,"X"]-150),
c(dat[i,"Y"]-150,dat[i,"Y"]+150,dat[i,"Y"]+150,dat[i,"Y"]-150,dat[i,"Y"]-150)))

Srs1 = Polygons(list(Sr1),paste("Number.",i,sep=""))

i=2
Sr2 = Polygon(cbind(c(dat[i,"X"]-150,dat[i,"X"] -150,dat[i,"X"]+150,dat[i,"X"]+150,dat[i,"X"]-150),
c(dat[i,"Y"]-150,dat[i,"Y"]+150,dat[i,"Y"]+150,dat[i,"Y"]-150,dat[i,"Y"]-150)))

Srs2 = Polygons(list(Sr2),paste("Number.",i,sep=""))

SpP = SpatialPolygons(list(Srs1, Srs2), 1:2)

SrDf = SpatialPolygonsDataFrame(SpP, dat[1:2,])

writePolyShape(SrDf, fn="TvaRutor")

However, my data frame "dat" has ~10000 coords and want to write a loop for doing what is done above. My
problem is how to store Sr1, Srs1, Sr2, Srs2,... Sr10000, Srs10000 etc to be used in SpatialPolygons().

Help on this is much appreciated.

```

1 Jul 2010 11:31

### Re: make thousands SpatialPolygons()

```Thanks Roman, this worked fine:
n=5
my.list <- vector('list', n)
for(i in 1:n) {
my.list[[i]] = Polygons(list(Polygon(cbind(c(dat[i,"X"]-150,dat[i,"X"]-150,dat[i,"X"]+150,dat[i,"X"]+150,dat[i,"X"]-150),
c(dat[i,"Y"]-150,dat[i,"Y"]+150,dat[i,"Y"]+150,dat[i,"Y"]-150,dat[i,"Y"]-150)))),paste("Number.",i,sep=""))
}
SpP = SpatialPolygons(my.list, 1:n)
SrDf = SpatialPolygonsDataFrame(SpP, dat[1:n,])
writePolyShape(SrDf, fn="FemRutor")

Tord

From: romunov [mailto:romunov <at> gmail.com]
Sent: den 1 juli 2010 11:14
To: Tord Snäll
Subject: Re: [R-sig-Geo] make thousands SpatialPolygons()

I would make a list out of all Sr/Srs objects - you will need it to make Polygons anyway.

I work with lists according to The R
Inferno<http://lib.stat.cmu.edu/S/Spoetry/Tutor/R_inferno.pdf> something along the lines:
my.list <- vector('list', n)
for(i in 1:n) {
my.list[[i]] <- your_object
}

Cheers,
Roman
```

1 Jul 2010 16:05

### Gaussian grids in R

```Dear  R-sig-geo members,

I wonder if any of you have suggestions about how to import and use
rasters with a spectral Gaussian T62 CRS in R.  I'm not very familiar
with this coordinate system (commonly used for climate models), and I
haven't been able to discover tools for dealing with them in R, either
in PROJ.4 or mapproj or anywhere else, despite extensive scouring of the
web.  Details of Gaussian grids can be seen at
http://en.wikipedia.org/wiki/Gaussian_grid and
http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID98

Specifically, I'm trying to use output from the NCEP CFS climate models
(see http://cfs.ncep.noaa.gov/cfs/monthly/tmp2m/ -- these are in GRIB
format).  A different example of a T62 grid can be downloaded in netCDF
format at http://www.sage.wisc.edu/iamdata/grid.php.

The goal is to transform the CRS into geographic lat/lon coordinates to
be more compatible with our other datasets.  The big problem with
getting these data into a SpatialGrid or SpatialPixelsDataFrame is that
the latitude cells are not equally spaced - they get smaller towards the
poles.  I may be able to approximate it to a regular grid by ignoring
the changes to the poles (they are rather small after all), but I'd
rather do it properly if possible.

Code below shows extracting latitudes from a T62 grid and plotting the
changes in cell size with latitude.

library('ncdf')
library('sp')
```

1 Jul 2010 16:53

### Nugget from variogram.fit vs GLS

```Hi All,

I'm verifying coefficients generated with GSTAT kriging with external
drift using GLS. The coefficients match absolutely perfectly. The range
also matches perfectly. But there is a big difference in what GLS gives
as the output for the nugget. Can anyone explain why this is?

Zev

###### KED

# here's where I fit the variogram and the fit output

v.fit = fit.variogram(v, vgm(1.5, "Exp", 20, 0.5))

> v.fit
model    psill     range
1   Nug 1.302310 0.0000000
2   Exp 1.413246 0.8098587

# the fit is used in a KED model
g = gstat(formula=formula(MODEL), data = so2, model = v.fit)

# here I get the coefs (these coefs match those from GLS below perfectly)
UKcoef<-predict(g, lattice, BLUE=c(TRUE,TRUE))

###### GLS

# use the same parameters as above and fix them

```

1 Jul 2010 17:40

### spTransform changes the class

```Hi!

I use spTransform to reproject:

> class(eugrd025)
[1] "SpatialGrid"
attr(,"package")
[1] "sp"

> projection(eugrd025)
[1] "+proj=longlat +ellps=WGS84"

> eugrd025EFDC <- spTransform(eugrd025, CRSobj=CRS(projection(eugrd025EFDC_SPDF)))

> projection(eugrd025EFDC)
[1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
+ellps=GRS80 +units=m +no_defs"

This is good, but:
> class(eugrd025EFDC)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"

Why is the class changed to SpatialPoints? I need a grid object for
> eugrd025EFDCr <- raster(eugrd025EFDC)
> Br025 <- resample(Br,eugrd025EFDCr)

don't I?

```

1 Jul 2010 17:46

### Re: spTransform changes the class

```spTransform cannot reproject a grid - this will (usually) requires
destructive resampling of the data to a completely new grid. There are
functions in the raster package to do this: projectRaster.

I think this would work, but there may be important details that you
will need to investigate:

raster.eugrd025 <- raster(eugrd025)

eugrd025EFDC <- projectRaster(raster.eugrd025,
CRSobj=CRS(projection(eugrd025EFDC_SPDF)))

On Fri, Jul 2, 2010 at 1:40 AM, Agustin Lobo <alobolistas <at> gmail.com> wrote:
> Hi!
>
> I use spTransform to reproject:
>
>> class(eugrd025)
> [1] "SpatialGrid"
> attr(,"package")
> [1] "sp"
>
>> projection(eugrd025)
> [1] "+proj=longlat +ellps=WGS84"
>
>> eugrd025EFDC <- spTransform(eugrd025, CRSobj=CRS(projection(eugrd025EFDC_SPDF)))
>
>> projection(eugrd025EFDC)
> [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
> +ellps=GRS80 +units=m +no_defs"
```

1 Jul 2010 17:48

### Re: spTransform changes the class

```Ah, sorry that was completely wrong. I think this is right, but still untested:

raster.eugrd025 <- raster(eugrd025)

pr <- projectExtent(raster.eugrd025, projection(eugrd025EFDC_SPDF))

eugrd025EFDC <- projectRaster(raster.eugrd025, pr)

On Fri, Jul 2, 2010 at 1:46 AM, Michael Sumner <mdsumner <at> gmail.com> wrote:
> spTransform cannot reproject a grid - this will (usually) requires
> destructive resampling of the data to a completely new grid. There are
> functions in the raster package to do this: projectRaster.
>
> I think this would work, but there may be important details that you
> will need to investigate:
>
> raster.eugrd025 <- raster(eugrd025)
>
> eugrd025EFDC <- projectRaster(raster.eugrd025,
> CRSobj=CRS(projection(eugrd025EFDC_SPDF)))
>
> On Fri, Jul 2, 2010 at 1:40 AM, Agustin Lobo <alobolistas <at> gmail.com> wrote:
>> Hi!
>>
>> I use spTransform to reproject:
>>
>>> class(eugrd025)
>> [1] "SpatialGrid"
>> attr(,"package")
>> [1] "sp"
```

1 Jul 2010 17:51

### Re: spTransform changes the class

```Hi!

Are you sure that the *spTransform *method can be applied to a SpatialGrid
object? I think it only works for SpatialPoints, SpatialLines and
SpatialPolygons.
If you are working with rasters, perhaps you should use the function
*projectRaster
*from the package raster.

I hope it helps...

Isa.

2010/7/1 Agustin Lobo <alobolistas <at> gmail.com>

> Hi!
>
> I use spTransform to reproject:
>
> > class(eugrd025)
> [1] "SpatialGrid"
> attr(,"package")
> [1] "sp"
>
> > projection(eugrd025)
> [1] "+proj=longlat +ellps=WGS84"
>
> > eugrd025EFDC <- spTransform(eugrd025,
> CRSobj=CRS(projection(eugrd025EFDC_SPDF)))
>
```

1 Jul 2010 18:00

### Re: spTransform changes the class

```I forgot to mention the first time, spTransform will not reproject a
grid but it will coerce the grid to points and then return reprojected
versions of those. That's why Agustin sees the class change from grid
to points.

You can see this in the method definitions for SpatialGrid/PixelDataFrame, e.g.

setMethod("spTransform", signature("SpatialPixelsDataFrame", "CRS"),
function(x, CRSobj, ...)
spTransform(as(x, "SpatialPointsDataFrame"), CRSobj, ...))

setMethod("spTransform", signature("SpatialGridDataFrame", "CRS"),
function(x, CRSobj, ...)
spTransform(as(x, "SpatialPixelsDataFrame"), CRSobj, ...))

Cheers, Mike.

On Fri, Jul 2, 2010 at 1:51 AM, isabelle boulangeat
<isabelle.boulangeat <at> gmail.com> wrote:
> Hi!
>
> Are you sure that the *spTransform *method can be applied to a SpatialGrid
> object? I think it only works for SpatialPoints, SpatialLines and
> SpatialPolygons.
> If you are working with rasters, perhaps you should use the function
> *projectRaster
> *from the package raster.
>
> I hope it helps...
>
```