headers
Agustin Lobo | 1 Jul 2010 07:02
Picon

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:
https://launchpad.net/~ubuntugis/+archive/ubuntugis-unstable
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
(Continue reading)

Permalink | Reply | 0 comments |
headers
Tord Snäll | 1 Jul 2010 10:59
Picon
Picon
Favicon

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.
(Continue reading)

Permalink | Reply | 0 comments |
headers
Tord Snäll | 1 Jul 2010 11:31
Picon
Picon
Favicon

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) {
... additional manipulation...
my.list[[i]] <- your_object
}

Cheers,
Roman
(Continue reading)

Permalink | Reply | 0 comments |
headers
Matthew Landis | 1 Jul 2010 16:05
Favicon

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')
nc <- 'C:/scratch/cfs/pop_den_1995_T62.nc' # Downloaded  from
(Continue reading)

Permalink | Reply | 5 comments |
headers
Zev Ross | 1 Jul 2010 16:53
Favicon

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
(Continue reading)

Permalink | Reply | 1 comment |
headers
Agustin Lobo | 1 Jul 2010 17:40
Picon

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?
(Continue reading)

Permalink | Reply | 3 comments |
headers
Michael Sumner | 1 Jul 2010 17:46
Picon
Gravatar

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"
(Continue reading)

Permalink | Reply | 0 comments |
headers
Michael Sumner | 1 Jul 2010 17:48
Picon
Gravatar

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"
(Continue reading)

Permalink | Reply | 3 comments |
headers
isabelle boulangeat | 1 Jul 2010 17:51
Picon

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)))
>
(Continue reading)

Permalink | Reply | 1 comment |
headers
Michael Sumner | 1 Jul 2010 18:00
Picon
Gravatar

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...
>
(Continue reading)

Permalink | Reply | 0 comments |

Gmane