[R-sig-Geo] A new ASTER Global DEM data set
Yong Li
yong.li at unimelb.edu.au
Fri Nov 20 12:39:25 CET 2009
Dear all folks,
I was informed in the 6th Digital Earth Conference that there is a better place to acquire high resolution of global DEM developed by ASTER, called GDEM with 30 m resolution, and fantastically free of charge. I tried some here (http://www.gdem.aster.ersdac.or.jp/) and it is really better than SRTM if you are outside USA.
Hope you will enjoy the free meal.
Cheers,
Yong
-----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 r-sig-geo-request at stat.math.ethz.ch
Sent: 2009年11月20日 22:00
To: r-sig-geo at stat.math.ethz.ch
Subject: R-sig-Geo Digest, Vol 75, Issue 18
Send R-sig-Geo mailing list submissions to
r-sig-geo at stat.math.ethz.ch
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
r-sig-geo-request at stat.math.ethz.ch
You can reach the person managing the list at
r-sig-geo-owner at stat.math.ethz.ch
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-Geo digest..."
Today's Topics:
1. Re: thinning a SpatialPolygonsDataFrame (Roger Bivand)
2. Spatial3dArray - coordinates method (Torleif Markussen Lunde)
3. Re: Spatial3dArray - coordinates method (Edzer Pebesma)
4. Re: thinning a SpatialPolygonsDataFrame (Michael Friendly)
5. Filling in holes in DTM (Tobin Cara)
6. Re: Filling in holes in DTM (milton ruser)
7. Re: Filling in holes in DTM (Tomislav Hengl)
8. execGRASS() for r.mapcalc? (Agustin Lobo)
9. Lambert projection to lat/long question.... (Blair Christian)
10. Re: Filling in holes in DTM (Robert J. Hijmans)
11. Re: Lambert projection to lat/long question....
(Robert J. Hijmans)
12. Re: Lambert projection to lat/long question....
(Robert J. Hijmans)
13. flags in execGRASS (Agustin Lobo)
14. Re: Lambert projection to lat/long question.... (Blair Christian)
15. Re: execGRASS() for r.mapcalc? (Roger Bivand)
16. Re: flags in execGRASS (Roger Bivand)
17. Re: Spatial3dArray - coordinates method (Michael Sumner)
18. Re: Spatial3dArray - coordinates method (Torleif Markussen Lunde)
19. Problem using mahasuhab (Ned Horning)
20. Re: Problem using mahasuhab (Mathieu Basille)
21. Re: Problem using mahasuhab (Ned Horning)
22. changing the data type of a gdal dataset (Alan Swanson)
----------------------------------------------------------------------
Message: 1
Date: Thu, 19 Nov 2009 13:02:11 +0100 (CET)
From: Roger Bivand <Roger.Bivand at nhh.no>
Subject: Re: [R-sig-Geo] thinning a SpatialPolygonsDataFrame
To: Pinaud David <pinaud at cebc.cnrs.fr>
Cc: r-sig-geo at stat.math.ethz.ch, Michael Friendly <friendly at yorku.ca>
Message-ID: <alpine.LRH.2.00.0911191253370.28000 at reclus.nhh.no>
Content-Type: text/plain; charset="iso-8859-1"; Format="flowed"
On Thu, 19 Nov 2009, Pinaud David wrote:
> Hi Michael,
> Roger is fully right, this function does not preserve the topology, so be
> aware of that, some problems can occur...
> If you want to use shapefiles::dp() just for raw plotting and visual
> simplification, you can try (on the fly):
>
> library(rgdal)
> library(shapefiles)
>
> fr <- readOGR("polygonFRA_WGS84.shp", "polygonFRA_WGS84") # a shapefile of
> France with complex topology (holes, islands...) in WGS84 coordinates
> pp <- slot(fr, "polygons") # take the polygons
> cf <- coordinates(slot(pp[[39]], "Polygons")[[1]]) # extract the coordinates
> of the main polygon ("continental France")
> pf <- list(x=cf[,1], y=cf[,2]) # list of coordinates, as dp() needs a list
> and not a matrix or dataframe...
> cf1 <- dp(pf, 0.1) # simplification, with a bandwith of 0.1 decimal degree
> plot(fr) # to see the result
> points(cf1, col="red", t="l")
This seems to work OK given that slivers are not important:
library(sp)
library(Guerry)
# installed from R-Forge
data(gfrance)
object.size(gfrance)
pls <- slot(gfrance, "polygons")
pls_dp <- vector(mode="list", length=length(pls))
require(shapefiles)
tol <- 2500
minarea <- 500000
for (i in 1:length(pls)) {
Pls <- slot(pls[[i]], "Polygons")
Pls_dp <- vector(mode="list", length=length(Pls))
for (j in 1:length(Pls)) {
crds <- slot(Pls[[j]], "coords")
crds_s <- dp(list(x=crds[,1], y=crds[,2]), tolerance=tol)
crds_s <- do.call("cbind", crds_s)
if(!identical(crds_s[1,], crds_s[nrow(crds_s),]))
crds_s <- rbind(crds_s, crds_s[1,])
Pls_dp[[j]] <- Polygon(crds_s)
}
Keep <- logical(length(Pls_dp))
for (j in 1:length(Pls_dp)) {
Keep[j] <- TRUE
if (slot(Pls_dp[[j]], "area") < minarea) Keep[j] <- FALSE
}
Pls_dp <- Pls_dp[Keep]
pls_dp[[i]] <- Polygons(Pls_dp, ID=slot(pls[[i]], "ID"))
}
gfrance_dp <- SpatialPolygonsDataFrame(SpatialPolygons(pls_dp),
data=slot(gfrance, "data"))
object.size(gfrance_dp)
If this was a function, the input would be an object inheriting from
SpatialPolygons, the tolerance, and a minimum polygon area. The removal of
small Polygon objects needs protecting from removing all belonging to a
given Polygons object. The CRS also needs copying across. The objects are
rebuilt to correct areas, centroids, plot orders, and the bounding box,
all of which may change. For figures, the companion thread on R-help is
relevant, PDF output for choropleth maps is often a good deal larger in
file size than that of the equivalent PNG device.
I'll ask the shapefiles authors whether I can copy dp() to maptools and
include suitable methods - this will help in benchmarking the future GEOS
version.
Hope this helps,
Roger
>
> HTH
> David
>
> Michael Friendly a ?crit :
>> I think, for my application, I'd be happy with the D-P polyline
>> simplification algorithm, because that is what is used
>> in SAS (worked well), and I don't think there are unusual topologies in my
>> map of France that would be so severely
>> out of whack as to lead to *significant* visual artifacts. In fact, you
>> might well expect some artifacts from any visual
>> thinning, but it's a matter of the tradeoff in the way the thinned map is
>> used in a visualization. Mark Monmonier's
>> US Visibility Map might be an extremely thinned, but highly useful example.
>>
>> For R spatial analysis, I think this is worth pursuing and integrating into
>> sp methods. In SAS, proc greduce works simply
>> by adding another variable, density, to the (x,y) coordinates of the
>> spatial polygons, density %in% 1:5, where density==5
>> is the full map. It is then a simple matter to subset the polygon outlines
>> by saying
>>
>> data smallmap;
>> set mymap;
>> where density<4;
>>
>> or
>>
>> proc gmap map=mymap(where=(density<4));
>> ...
>>
>> Meanwhile, I can't see easily how I could use shapefiles::dp() to thin my
>> Guerry::gfrance maps, because the documentation is,
>> shall we say, somewhat thin.
>> -Michael
>>
>>
>>
>> Roger Bivand wrote:
>>> On Wed, 18 Nov 2009, Pinaud David wrote:
>>>
>>>> Hi Michael,
>>>> maybe you should try the function dp() in the package shapefiles that is
>>>> an implementation of the Douglas-Peucker polyLine simplification
>>>> algorithm.
>>>
>>> Note that its help page does warn that it is not topology-preserving, that
>>> is that lines are generalised, but that coincident lines (boundaries of
>>> neighbouring polygons) may be generalised differently. GEOS offers a
>>> topology-preserving line generalisation facility, which ought to take
>>> longer but do better than dp(), because it will not lead to visual
>>> artefacts (overlapping polygons, interpolygon slivers, etc.).
>>>
>>> Roger
>>>
>>>> HTH
>>>> David
>>>>
>>>> Michael Friendly a ?crit :
>>>>> The Guerry package contains two maps of france (gfrance, gfrance85)
>>>>> which are quite detailed and large in size (‾900K).
>>>>> In writing a vignette for the package, there are quite a few figures
>>>>> that use the map multiple times in a layout, and
>>>>> consequently result in huge file sizes for the .PDF files created. For
>>>>> these purposes, the map need not be nearly
>>>>> so detailed.
>>>>>
>>>>> I'm wondering if there is a facility to "thin" the map by drawing it at
>>>>> a lower density of lines in the polygon regions.
>>>>> When I was working with SAS, there was a GREDUCE procedure that did this
>>>>> nicely.
>>>>>
>>>>> thanks,
>>>>> -Michael
>>>>>
>>>>
>>>>
>>>
>>
>>
>
>
--
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
------------------------------
Message: 2
Date: Thu, 19 Nov 2009 13:28:55 +0100
From: Torleif Markussen Lunde <torleif.lunde at cih.uib.no>
Subject: [R-sig-Geo] Spatial3dArray - coordinates method
To: r-sig-geo at stat.math.ethz.ch
Message-ID: <200911191328.55628.torleif.lunde at cih.uib.no>
Content-Type: Text/Plain; charset="iso-8859-15"
Hi
To read netcdf data (or any other "gridded" spatial time data) I find it
convenient to define new classes Spatial3dArray and Spatial4dArray.
setClass("Spatial3dArray",
representation("Spatial", data = "array", coords = "list",
time = "character", btime = "character"),
prototype= list(data = array(NA, c(1,1,1,1)),
bbox=matrix(NA),
proj4string = CRS(as.character(NA)),
coords = list(1,1),
time = "posix",
btime = "posix"))
##########################################
###################EXAMPLE##################
##########################################
x <- matrix(seq(-10, 10, length = 100), 100, 100,
byrow = FALSE)
y <- matrix(seq(-10, 10, length = 100), 100, 100,
byrow = TRUE)
tm <- 1:10
tm.c <- as.character(seq(as.POSIXct("2002-01-01 06:00:00",
"2002-01-01 15:00:00"),
by="hours",
length.out=10))
z <- array(NA, c(dim(x)[1], dim(x)[2], length(tm.c), 1))
for (i in 1:10) {
z[,,i,] <- i * ( sin(sqrt(x^2+y^2)))
}
sin3dA <- new("Spatial3dArray",
data = z,
coords = list(x, y),
bbox = matrix(c(min(x), min(y), max(x), max(y), 2, 2), 2, 2,
dimnames = list(NULL, c("min","max"))),
time = tm.c,
btime = c(min(tm.c), max(tm.c)))
dimnames(slot(sin3dA, "data")) = list(NULL,
NULL,
slot(sin3dA, "time"),
c("a"))
names(slot(sin3dA, "coords")) <- c("x", "y")
##########################################
for the coordinates method I would like to have two options on how to return
the coordinates; "list" or default "sp":
coordinates.3dArray <- function (obj, type = "sp") {
lat <- slot(obj, "coords")[[1]]
long <- slot(obj, "coords")[[2]]
if (type == "list") {
return(list(long=long, lat=lat))
} else if (type == "sp") {
res <- as.matrix(cbind(c(long), c(lat)))
dimnames(res) <- list(NULL, c("x1", "x2"))
}
}
setMethod("coordinates", signature("Spatial3dArray"), coordinates.3dArray)
This means that the default coordinates method in sp has to include the option
"..." . Would it be possible to include this in a future release of sp?
The reason I want to keep the list option is to use a matrix oriented approach
in spplot, overlay, etc. methods. I also feel having a matrix/array approach
with these kind of data makes sense. Allowing type = "sp" means overlay() will
work more or less out of the box (however I would like to return a matrix),
and still I could get the list/matrix when desired.
Best wishes
Torleif Markussen Lunde
Centre for International Health
Bjerknes Centre for Climate Research
University of Bergen
Norway
------------------------------
Message: 3
Date: Thu, 19 Nov 2009 14:26:35 +0100
From: Edzer Pebesma <edzer.pebesma at uni-muenster.de>
Subject: Re: [R-sig-Geo] Spatial3dArray - coordinates method
To: Torleif Markussen Lunde <torleif.lunde at cih.uib.no>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4B05478B.3080800 at uni-muenster.de>
Content-Type: text/plain; charset=ISO-8859-1
No problems; sp in csv now has this, the next release will have it.
Torleif Markussen Lunde wrote:
> Hi
>
> To read netcdf data (or any other "gridded" spatial time data) I find it
> convenient to define new classes Spatial3dArray and Spatial4dArray.
>
> setClass("Spatial3dArray",
> representation("Spatial", data = "array", coords = "list",
> time = "character", btime = "character"),
> prototype= list(data = array(NA, c(1,1,1,1)),
> bbox=matrix(NA),
> proj4string = CRS(as.character(NA)),
> coords = list(1,1),
> time = "posix",
> btime = "posix"))
>
>
> ##########################################
> ###################EXAMPLE##################
> ##########################################
>
> x <- matrix(seq(-10, 10, length = 100), 100, 100,
> byrow = FALSE)
> y <- matrix(seq(-10, 10, length = 100), 100, 100,
> byrow = TRUE)
>
> tm <- 1:10
> tm.c <- as.character(seq(as.POSIXct("2002-01-01 06:00:00",
> "2002-01-01 15:00:00"),
> by="hours",
> length.out=10))
>
> z <- array(NA, c(dim(x)[1], dim(x)[2], length(tm.c), 1))
>
> for (i in 1:10) {
> z[,,i,] <- i * ( sin(sqrt(x^2+y^2)))
> }
>
> sin3dA <- new("Spatial3dArray",
> data = z,
> coords = list(x, y),
> bbox = matrix(c(min(x), min(y), max(x), max(y), 2, 2), 2, 2,
> dimnames = list(NULL, c("min","max"))),
> time = tm.c,
> btime = c(min(tm.c), max(tm.c)))
>
> dimnames(slot(sin3dA, "data")) = list(NULL,
> NULL,
> slot(sin3dA, "time"),
> c("a"))
> names(slot(sin3dA, "coords")) <- c("x", "y")
>
>
> ##########################################
>
> for the coordinates method I would like to have two options on how to return
> the coordinates; "list" or default "sp":
>
> coordinates.3dArray <- function (obj, type = "sp") {
> lat <- slot(obj, "coords")[[1]]
> long <- slot(obj, "coords")[[2]]
> if (type == "list") {
> return(list(long=long, lat=lat))
> } else if (type == "sp") {
> res <- as.matrix(cbind(c(long), c(lat)))
> dimnames(res) <- list(NULL, c("x1", "x2"))
> }
> }
> setMethod("coordinates", signature("Spatial3dArray"), coordinates.3dArray)
>
> This means that the default coordinates method in sp has to include the option
> "..." . Would it be possible to include this in a future release of sp?
>
> The reason I want to keep the list option is to use a matrix oriented approach
> in spplot, overlay, etc. methods. I also feel having a matrix/array approach
> with these kind of data makes sense. Allowing type = "sp" means overlay() will
> work more or less out of the box (however I would like to return a matrix),
> and still I could get the list/matrix when desired.
>
>
> Best wishes
> Torleif Markussen Lunde
> Centre for International Health
> Bjerknes Centre for Climate Research
> University of Bergen
> Norway
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of M?nster
Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
------------------------------
Message: 4
Date: Thu, 19 Nov 2009 08:55:10 -0500
From: Michael Friendly <friendly at yorku.ca>
Subject: Re: [R-sig-Geo] thinning a SpatialPolygonsDataFrame
To: Roger.Bivand at nhh.no
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4B054E3E.4060005 at yorku.ca>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Thanks very much for doing the heavy lifting here, Roger.
For use in the Guerry package, I've turned this into a function (rather
than saving other versions of gfrance).
If you include something like this in sp, I'll withdraw it from Guerry.
thinnedSpatialPoly <- function(SP, tolerance, minarea) {
if (!require(shapefiles)) stop("shapefiles package is required")
stopifnot(inherits(SP, "SpatialPolygons"))
# TODO: determine and set defaults for tolerance, minarea
# TODO: suppress warnings: "In Polygon(crds_s) : Non-finite label
point detected and replaced"
pls <- slot(SP, "polygons")
pls_dp <- vector(mode="list", length=length(pls))
for (i in 1:length(pls)) {
Pls <- slot(pls[[i]], "Polygons")
Pls_dp <- vector(mode="list", length=length(Pls))
for (j in 1:length(Pls)) {
crds <- slot(Pls[[j]], "coords")
crds_s <- dp(list(x=crds[,1], y=crds[,2]), tolerance=tolerance)
crds_s <- do.call("cbind", crds_s)
if(!identical(crds_s[1,], crds_s[nrow(crds_s),]))
crds_s <- rbind(crds_s, crds_s[1,])
Pls_dp[[j]] <- Polygon(crds_s)
}
Keep <- logical(length(Pls_dp))
for (j in 1:length(Pls_dp)) {
Keep[j] <- TRUE
if (slot(Pls_dp[[j]], "area") < minarea) Keep[j] <- FALSE
}
Pls_dp <- Pls_dp[Keep]
pls_dp[[i]] <- Polygons(Pls_dp, ID=slot(pls[[i]], "ID"))
}
SP_dp <- SpatialPolygons(pls_dp, proj4string=slot(SP, "proj4string"))
if(inherits(SP, "SpatialPolygonsDataFrame")) {
data <- slot(SP, "data")
SP_dp <- SpatialPolygonsDataFrame(SP_dp, data=data)
}
SP_dp
}
best,
-Michael
Roger Bivand wrote:
> On Thu, 19 Nov 2009, Pinaud David wrote:
>
>> Hi Michael,
>> Roger is fully right, this function does not preserve the topology,
>> so be aware of that, some problems can occur...
>> If you want to use shapefiles::dp() just for raw plotting and visual
>> simplification, you can try (on the fly):
>>
>> library(rgdal)
>> library(shapefiles)
>>
>> fr <- readOGR("polygonFRA_WGS84.shp", "polygonFRA_WGS84") # a
>> shapefile of France with complex topology (holes, islands...) in
>> WGS84 coordinates
>> pp <- slot(fr, "polygons") # take the polygons
>> cf <- coordinates(slot(pp[[39]], "Polygons")[[1]]) # extract the
>> coordinates of the main polygon ("continental France")
>> pf <- list(x=cf[,1], y=cf[,2]) # list of coordinates, as dp() needs a
>> list and not a matrix or dataframe...
>> cf1 <- dp(pf, 0.1) # simplification, with a bandwith of 0.1 decimal
>> degree
>> plot(fr) # to see the result
>> points(cf1, col="red", t="l")
>
> This seems to work OK given that slivers are not important:
> library(sp)
> library(Guerry)
> # installed from R-Forge
> data(gfrance)
> object.size(gfrance)
> pls <- slot(gfrance, "polygons")
> pls_dp <- vector(mode="list", length=length(pls))
> require(shapefiles)
> tol <- 2500
> minarea <- 500000
> for (i in 1:length(pls)) {
> Pls <- slot(pls[[i]], "Polygons")
> Pls_dp <- vector(mode="list", length=length(Pls))
> for (j in 1:length(Pls)) {
> crds <- slot(Pls[[j]], "coords")
> crds_s <- dp(list(x=crds[,1], y=crds[,2]), tolerance=tol)
> crds_s <- do.call("cbind", crds_s)
> if(!identical(crds_s[1,], crds_s[nrow(crds_s),]))
> crds_s <- rbind(crds_s, crds_s[1,])
> Pls_dp[[j]] <- Polygon(crds_s)
> }
> Keep <- logical(length(Pls_dp))
> for (j in 1:length(Pls_dp)) {
> Keep[j] <- TRUE
> if (slot(Pls_dp[[j]], "area") < minarea) Keep[j] <- FALSE
> }
> Pls_dp <- Pls_dp[Keep]
> pls_dp[[i]] <- Polygons(Pls_dp, ID=slot(pls[[i]], "ID"))
> }
> gfrance_dp <- SpatialPolygonsDataFrame(SpatialPolygons(pls_dp),
> data=slot(gfrance, "data"))
> object.size(gfrance_dp)
>
> If this was a function, the input would be an object inheriting from
> SpatialPolygons, the tolerance, and a minimum polygon area. The
> removal of small Polygon objects needs protecting from removing all
> belonging to a given Polygons object. The CRS also needs copying
> across. The objects are rebuilt to correct areas, centroids, plot
> orders, and the bounding box, all of which may change. For figures,
> the companion thread on R-help is relevant, PDF output for choropleth
> maps is often a good deal larger in file size than that of the
> equivalent PNG device.
>
> I'll ask the shapefiles authors whether I can copy dp() to maptools
> and include suitable methods - this will help in benchmarking the
> future GEOS version.
>
> Hope this helps,
>
> Roger
>
>>
>> HTH
>> David
>>
>> Michael Friendly a ?crit :
>>> I think, for my application, I'd be happy with the D-P polyline
>>> simplification algorithm, because that is what is used
>>> in SAS (worked well), and I don't think there are unusual topologies
>>> in my map of France that would be so severely
>>> out of whack as to lead to *significant* visual artifacts. In fact,
>>> you might well expect some artifacts from any visual
>>> thinning, but it's a matter of the tradeoff in the way the thinned
>>> map is used in a visualization. Mark Monmonier's
>>> US Visibility Map might be an extremely thinned, but highly useful
>>> example.
>>>
>>> For R spatial analysis, I think this is worth pursuing and
>>> integrating into sp methods. In SAS, proc greduce works simply
>>> by adding another variable, density, to the (x,y) coordinates of the
>>> spatial polygons, density %in% 1:5, where density==5
>>> is the full map. It is then a simple matter to subset the polygon
>>> outlines by saying
>>>
>>> data smallmap;
>>> set mymap;
>>> where density<4;
>>>
>>> or
>>>
>>> proc gmap map=mymap(where=(density<4));
>>> ...
>>>
>>> Meanwhile, I can't see easily how I could use shapefiles::dp() to
>>> thin my Guerry::gfrance maps, because the documentation is,
>>> shall we say, somewhat thin.
>>> -Michael
>>>
>>>
>>>
>>> Roger Bivand wrote:
>>>> On Wed, 18 Nov 2009, Pinaud David wrote:
>>>>
>>>>> Hi Michael,
>>>>> maybe you should try the function dp() in the package shapefiles
>>>>> that is an implementation of the Douglas-Peucker polyLine
>>>>> simplification algorithm.
>>>>
>>>> Note that its help page does warn that it is not
>>>> topology-preserving, that is that lines are generalised, but that
>>>> coincident lines (boundaries of neighbouring polygons) may be
>>>> generalised differently. GEOS offers a topology-preserving line
>>>> generalisation facility, which ought to take longer but do better
>>>> than dp(), because it will not lead to visual artefacts
>>>> (overlapping polygons, interpolygon slivers, etc.).
>>>>
>>>> Roger
>>>>
>>>>> HTH
>>>>> David
>>>>>
>>>>> Michael Friendly a ?crit :
>>>>>> The Guerry package contains two maps of france (gfrance,
>>>>>> gfrance85) which are quite detailed and large in size (‾900K).
>>>>>> In writing a vignette for the package, there are quite a few
>>>>>> figures that use the map multiple times in a layout, and
>>>>>> consequently result in huge file sizes for the .PDF files
>>>>>> created. For these purposes, the map need not be nearly
>>>>>> so detailed.
>>>>>>
>>>>>> I'm wondering if there is a facility to "thin" the map by drawing
>>>>>> it at a lower density of lines in the polygon regions.
>>>>>> When I was working with SAS, there was a GREDUCE procedure that
>>>>>> did this nicely.
>>>>>>
>>>>>> thanks,
>>>>>> -Michael
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>
>>
>
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT M3J 1P3 CANADA
------------------------------
Message: 5
Date: Thu, 19 Nov 2009 15:16:23 +0100
From: Tobin Cara <cara.tobin at epfl.ch>
Subject: [R-sig-Geo] Filling in holes in DTM
To: "r-sig-geo at stat.math.ethz.ch" <r-sig-geo at stat.math.ethz.ch>
Message-ID:
<EDB94195DABE64488928DD39E53B8FC68BDFD3E8D9 at REX2.intranet.epfl.ch>
Content-Type: text/plain
Hello,
Would anyone know a good way to fill in holes within a DTM? There is no data inside these holes and it is affecting my calculations, so I prefer for an interpolation of the nearest neighbors to fill in a value or something similar.
Thank you,
Cara
[[alternative HTML version deleted]]
------------------------------
Message: 6
Date: Thu, 19 Nov 2009 10:04:58 -0500
From: milton ruser <milton.ruser at gmail.com>
Subject: Re: [R-sig-Geo] Filling in holes in DTM
To: Tobin Cara <cara.tobin at epfl.ch>
Cc: "r-sig-geo at stat.math.ethz.ch" <r-sig-geo at stat.math.ethz.ch>
Message-ID:
<3aaf1a030911190704s78a9fa2ah2fb3ee444be9a1e3 at mail.gmail.com>
Content-Type: text/plain
Hi Ton
Check r.watershed/r.terraflow/r.flow/r.fill modules on Grass.
cheers
miltinho
On Thu, Nov 19, 2009 at 9:16 AM, Tobin Cara <cara.tobin at epfl.ch> wrote:
> Hello,
>
> Would anyone know a good way to fill in holes within a DTM? There is no
> data inside these holes and it is affecting my calculations, so I prefer for
> an interpolation of the nearest neighbors to fill in a value or something
> similar.
>
> Thank you,
>
> Cara
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
[[alternative HTML version deleted]]
------------------------------
Message: 7
Date: Thu, 19 Nov 2009 16:14:56 +0100
From: "Tomislav Hengl" <hengl at spatial-analyst.net>
Subject: Re: [R-sig-Geo] Filling in holes in DTM
To: "'r-sig-geo'" <r-sig-geo at stat.math.ethz.ch>
Message-ID: <4117908DD26B4334AC97EE4600F89019 at pcibed193>
Content-Type: text/plain; charset="windows-1250"
Hi Cara,
There is a very efficient function in SAGA called "Close Gaps" that does exactly that. What makes it
especially efficient is that it allows you to set a mask map. See:
> rsaga.get.usage("grid_tools", 7)
SAGA CMD 2.0.3
library path: C:/Progra‾1/saga_vc/modules
library name: grid_tools
module name : Close Gaps
Usage: 7 -INPUT <str> [-MASK <str>] [-RESULT <str>] [-THRESHOLD <str>]
-INPUT:<str> Grid
Grid (input)
-MASK:<str> Mask
Grid (optional input)
-RESULT:<str> Changed Grid
Grid (optional output)
-THRESHOLD:<str> Tension Threshold
BR,
T. Hengl
http://home.medewerker.uva.nl/t.hengl/
> -----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 Tobin Cara
> Sent: Thursday, November 19, 2009 3:16 PM
> To: r-sig-geo at stat.math.ethz.ch
> Subject: [R-sig-Geo] Filling in holes in DTM
>
> Hello,
>
> Would anyone know a good way to fill in holes within a DTM? There is no data inside these holes
> and it is affecting my calculations, so I prefer for an interpolation of the nearest neighbors to
> fill in a value or something similar.
>
> Thank you,
>
> Cara
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
------------------------------
Message: 8
Date: Thu, 19 Nov 2009 18:09:57 +0100
From: Agustin Lobo <alobolistas at gmail.com>
Subject: [R-sig-Geo] execGRASS() for r.mapcalc?
To: sig-geo <r-sig-geo at stat.math.ethz.ch>
Message-ID: <4B057BE5.8040908 at gmail.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
I'm used to do:
> cmd = "r.mapcalc 'delme=100'"
> system(cmd)
and works, but how could I do the equivalent using execGRASS() ?
I've tried:
> doGRASS(cmd)
Error in parseGRASS(cmd) : r.mapcalc 'delme=100' not parsed
> execGRASS(cmd)
Error in parseGRASS(cmd) : r.mapcalc 'delme=100' not parsed
I can always use paste() to make any cmd string, but seems to me that
execGRASS() would be cleaner.
Thanks
Agus
------------------------------
Message: 9
Date: Thu, 19 Nov 2009 12:39:15 -0500
From: Blair Christian <blair.christian at gmail.com>
Subject: [R-sig-Geo] Lambert projection to lat/long question....
To: r-sig-geo at stat.math.ethz.ch
Message-ID:
<6c35a4fc0911190939o1ae8a7eeua800dcaeb4d3f21d at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Hi All,
Before trying to reinvent the wheel, I was wondering if the mapproj or
any other library had code to take lambert coordinates (along with the
parameters used to generate them) and convert them to lat long.? I
have grid data in lambert coords, and want to generate the long/lat
grid corners from them.? I noticed mapproj goes long/lat ->
projection, but couldn't find anything to go the other way. I was
going to use the reference in wolfram mathworld to do the inverse
transform, if anybody has any issues with their version.
http://mathworld.wolfram.com/LambertConformalConicProjection.html
Basically, I want to keep the data stored as grids, but every now and
then I want to create distance matrices or plot them (which I
currently do with an irregular grid class I made up), but don't want
to store polygons all the time when I can just store the basic grid.
All comments welcome.
Thanks,
Blair
------------------------------
Message: 10
Date: Thu, 19 Nov 2009 09:40:44 -0800
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
Subject: Re: [R-sig-Geo] Filling in holes in DTM
To: Tobin Cara <cara.tobin at epfl.ch>
Cc: "r-sig-geo at stat.math.ethz.ch" <r-sig-geo at stat.math.ethz.ch>
Message-ID:
<dc22b2570911190940u62269dbbo2392da32b0019a9d at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Tobin,
In the raster package you can use the focalNA function. It sets the
focal (neighborhood) value, according to a specified function (e.g.
mean), to cells with NA.
raster is available on R-Forge: install.packages("raster",
repos="http://R-Forge.R-project.org")
Robert
On Thu, Nov 19, 2009 at 6:16 AM, Tobin Cara <cara.tobin at epfl.ch> wrote:
> Hello,
>
> Would anyone know a good way to fill in holes within a DTM? There is no data inside these holes and it is affecting my calculations, so I prefer for an interpolation of the nearest neighbors to fill in a value or something similar.
>
> Thank you,
>
> Cara
>
> ? ? ? ?[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
------------------------------
Message: 11
Date: Thu, 19 Nov 2009 09:44:32 -0800
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
Subject: Re: [R-sig-Geo] Lambert projection to lat/long question....
To: Blair Christian <blair.christian at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID:
<dc22b2570911190944l49cafe4fx7b4367bf7e466e72 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Blair: Have a look at spTransform in rgdal. Robert
On Thu, Nov 19, 2009 at 9:39 AM, Blair Christian
<blair.christian at gmail.com> wrote:
> Hi All,
>
> Before trying to reinvent the wheel, I was wondering if the mapproj or
> any other library had code to take lambert coordinates (along with the
> parameters used to generate them) and convert them to lat long.? I
> have grid data in lambert coords, and want to generate the long/lat
> grid corners from them.? I noticed mapproj goes long/lat ->
> projection, but couldn't find anything to go the other way. ?I was
> going to use the reference in wolfram mathworld to do the inverse
> transform, if anybody has any issues with their version.
>
> http://mathworld.wolfram.com/LambertConformalConicProjection.html
>
> Basically, I want to keep the data stored as grids, but every now and
> then I want to create distance matrices or plot them (which I
> currently do with an irregular grid class I made up), but don't want
> to store polygons all the time when I can just store the basic grid.
> All comments welcome.
>
> Thanks,
> Blair
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
------------------------------
Message: 12
Date: Thu, 19 Nov 2009 09:46:13 -0800
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
Subject: Re: [R-sig-Geo] Lambert projection to lat/long question....
To: Blair Christian <blair.christian at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID:
<dc22b2570911190946w7f5d668bpa233ecb2763f24ad at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Blair, sorry, too quick:
Have a look at spTransform in rgdal and at projectRaster in the raster package.
Robert
On Thu, Nov 19, 2009 at 9:39 AM, Blair Christian
<blair.christian at gmail.com> wrote:
> Hi All,
>
> Before trying to reinvent the wheel, I was wondering if the mapproj or
> any other library had code to take lambert coordinates (along with the
> parameters used to generate them) and convert them to lat long.? I
> have grid data in lambert coords, and want to generate the long/lat
> grid corners from them.? I noticed mapproj goes long/lat ->
> projection, but couldn't find anything to go the other way. ?I was
> going to use the reference in wolfram mathworld to do the inverse
> transform, if anybody has any issues with their version.
>
> http://mathworld.wolfram.com/LambertConformalConicProjection.html
>
> Basically, I want to keep the data stored as grids, but every now and
> then I want to create distance matrices or plot them (which I
> currently do with an irregular grid class I made up), but don't want
> to store polygons all the time when I can just store the basic grid.
> All comments welcome.
>
> Thanks,
> Blair
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
------------------------------
Message: 13
Date: Thu, 19 Nov 2009 18:46:24 +0100
From: Agustin Lobo <alobolistas at gmail.com>
Subject: [R-sig-Geo] flags in execGRASS
To: sig-geo <r-sig-geo at stat.math.ethz.ch>
Message-ID: <4B058470.9000300 at gmail.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi!
The following works:
>
doGRASS("r.grow.distance",flags=c("overwrite"),parameters=list(distance="dummyd",input="dummy"))
[1] "r.grow.distance --overwrite distance=dummyd input=dummy"
But using "o" instead of "overwrite" does not
>
doGRASS("r.grow.distance",flags=c("o"),parameters=list(distance="dummyd",input="dummy"))
Error in doGRASS("r.grow.distance", flags = c("o"), parameters = list(distance =
"dummyd", :
Invalid flag value: o
despite the fact that the grass manual states:
> r.grow.distance --help
...(stuff deleted)
Usage:
r.grow.distance input=name [distance=name] [value=name]
[metric=string] [--overwrite] [--verbose] [--quiet]
Flags:
--o Allow output files to overwrite existing files
--v Verbose module output
--q Quiet module output
...(stuff deleted)
Could it be possible that execGRASS accept both forms for consistency with the
grass command? If this is not possible or appropriate, could a note be added to
the execGRASS manual page?
Thanks
Agus
------------------------------
Message: 14
Date: Thu, 19 Nov 2009 12:56:35 -0500
From: Blair Christian <blair.christian at gmail.com>
Subject: Re: [R-sig-Geo] Lambert projection to lat/long question....
To: "Robert J. Hijmans" <r.hijmans at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID:
<6c35a4fc0911190956x16b99092w2908ae9e0982765b at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Those are perfect. I wish I had the time to contribute as much
software as you do, Robert. Very impressive contributions on the
r-forge pages.
Thanks,
Blair
On Thu, Nov 19, 2009 at 12:46 PM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> Blair, sorry, too quick:
>
> Have a look at spTransform in rgdal and at projectRaster in the raster package.
>
> Robert
>
>
> On Thu, Nov 19, 2009 at 9:39 AM, Blair Christian
> <blair.christian at gmail.com> wrote:
>> Hi All,
>>
>> Before trying to reinvent the wheel, I was wondering if the mapproj or
>> any other library had code to take lambert coordinates (along with the
>> parameters used to generate them) and convert them to lat long.? I
>> have grid data in lambert coords, and want to generate the long/lat
>> grid corners from them.? I noticed mapproj goes long/lat ->
>> projection, but couldn't find anything to go the other way. ?I was
>> going to use the reference in wolfram mathworld to do the inverse
>> transform, if anybody has any issues with their version.
>>
>> http://mathworld.wolfram.com/LambertConformalConicProjection.html
>>
>> Basically, I want to keep the data stored as grids, but every now and
>> then I want to create distance matrices or plot them (which I
>> currently do with an irregular grid class I made up), but don't want
>> to store polygons all the time when I can just store the basic grid.
>> All comments welcome.
>>
>> Thanks,
>> Blair
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
------------------------------
Message: 15
Date: Thu, 19 Nov 2009 19:26:55 +0100 (CET)
From: Roger Bivand <Roger.Bivand at nhh.no>
Subject: Re: [R-sig-Geo] execGRASS() for r.mapcalc?
To: Agustin.Lobo at ija.csic.es
Cc: sig-geo <r-sig-geo at stat.math.ethz.ch>
Message-ID: <alpine.LRH.2.00.0911191921260.28832 at reclus.nhh.no>
Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII
On Thu, 19 Nov 2009, Agustin Lobo wrote:
> I'm used to do:
>> cmd = "r.mapcalc 'delme=100'"
>> system(cmd)
>
> and works, but how could I do the equivalent using execGRASS() ?
You cannot. r.mapcalc does not provide an --interface-description XML
self-definition. You have either to use system() as before (with different
quoting reimes for different OS), or r.mapcalculator in execGRASS(). This
problem is shared with other interfaces using --interface-description.
Roger
> I've tried:
>
>> doGRASS(cmd)
> Error in parseGRASS(cmd) : r.mapcalc 'delme=100' not parsed
>
>> execGRASS(cmd)
> Error in parseGRASS(cmd) : r.mapcalc 'delme=100' not parsed
>
> I can always use paste() to make any cmd string, but seems to me that
> execGRASS() would be cleaner.
>
> Thanks
>
> Agus
>
> _______________________________________________
> 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
------------------------------
Message: 16
Date: Thu, 19 Nov 2009 19:34:50 +0100 (CET)
From: Roger Bivand <Roger.Bivand at nhh.no>
Subject: Re: [R-sig-Geo] flags in execGRASS
To: Agustin.Lobo at ija.csic.es
Cc: sig-geo <r-sig-geo at stat.math.ethz.ch>
Message-ID: <alpine.LRH.2.00.0911191927160.28832 at reclus.nhh.no>
Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII
On Thu, 19 Nov 2009, Agustin Lobo wrote:
> Hi!
> The following works:
>>
> doGRASS("r.grow.distance",flags=c("overwrite"),parameters=list(distance="dummyd",input="dummy"))
> [1] "r.grow.distance --overwrite distance=dummyd input=dummy"
>
> But using "o" instead of "overwrite" does not
Because only "overwrite" is included in
r.grow.distance --interface-description
and that is what the flags are checked against. Unfortunately, many XML
descriptions are more demanding than the help pages, but for parseGRASS()
to work, it must use the self-definition without exception. execGRASS() is
designed for scripting, in which case keeping to tighter limits than the
help pages is OK. Some commands have both "o" and "overwrite" for
different things too, so confusion would ensue.
Roger
>>
> doGRASS("r.grow.distance",flags=c("o"),parameters=list(distance="dummyd",input="dummy"))
> Error in doGRASS("r.grow.distance", flags = c("o"), parameters =
> list(distance = "dummyd", :
> Invalid flag value: o
>
> despite the fact that the grass manual states:
>> r.grow.distance --help
> ...(stuff deleted)
>
> Usage:
> r.grow.distance input=name [distance=name] [value=name]
> [metric=string] [--overwrite] [--verbose] [--quiet]
>
> Flags:
> --o Allow output files to overwrite existing files
> --v Verbose module output
> --q Quiet module output
> ...(stuff deleted)
>
> Could it be possible that execGRASS accept both forms for consistency with
> the grass command? If this is not possible or appropriate, could a note be
> added to the execGRASS manual page?
>
> Thanks
>
> Agus
>
> _______________________________________________
> 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
------------------------------
Message: 17
Date: Fri, 20 Nov 2009 06:23:35 +1100
From: Michael Sumner <mdsumner at gmail.com>
Subject: Re: [R-sig-Geo] Spatial3dArray - coordinates method
To: Torleif Markussen Lunde <torleif.lunde at cih.uib.no>,
r-sig-geo at stat.math.ethz.ch
Message-ID:
<522664f80911191123u644807ccx7973ce8f734a0afe at mail.gmail.com>
Content-Type: text/plain; charset=UTF-8
Wow, this is great - I was thinking about this just yesterday.
Torleif: do you have an opinion on which NetCDF path is the most
useful for R with sp? RNetCDF or ncdf? GDAL is workable but takes
extra effort to build and then reconstruct 3d/4d from 2d bands. (I use
Windows mostly)
I use the RNetCDF package a lot, mainly because that's the one I
learnt to use first - there are binaries for Windows. It has some
problems in terms of R-style but they could be easily fixed.
Regards, Mike.
On Fri, Nov 20, 2009 at 12:26 AM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
> No problems; sp in csv now has this, the next release will have it.
>
> Torleif Markussen Lunde wrote:
>> Hi
>>
>> To read netcdf data (or any other "gridded" spatial time data) I find it
>> convenient to define new classes Spatial3dArray and Spatial4dArray.
>>
>> ?setClass("Spatial3dArray",
>> ? ? ? ? representation("Spatial", data = "array", coords = "list",
>> ? ? ? ? ? ? ? ? ? ? ? ? time = "character", btime = "character"),
>> ? ? ? ? prototype= list(data = array(NA, c(1,1,1,1)),
>> ? ? ? ? ? ? ? ? ? ? ? ? bbox=matrix(NA),
>> ? ? ? ? ? ? ? ? ? ? ? ? proj4string = CRS(as.character(NA)),
>> ? ? ? ? ? ? ? ? ? ? ? ? coords = list(1,1),
>> ? ? ? ? ? ? ? ? ? ? ? ? time = "posix",
>> ? ? ? ? ? ? ? ? ? ? ? ? btime = "posix"))
>>
>>
>> ##########################################
>> ###################EXAMPLE##################
>> ##########################################
>>
>> x <- matrix(seq(-10, 10, length = 100), 100, 100,
>> ? ? ? ? ? byrow = FALSE)
>> y <- matrix(seq(-10, 10, length = 100), 100, 100,
>> ? ? ? ? ? byrow = TRUE)
>>
>> tm <- 1:10
>> tm.c <- as.character(seq(as.POSIXct("2002-01-01 06:00:00",
>> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? "2002-01-01 15:00:00"),
>> ? ? ? ? ? ? ? ? ? ? ? ? by="hours",
>> ? ? ? ? ? ? ? ? ? ? ? ? length.out=10))
>>
>> z <- array(NA, c(dim(x)[1], dim(x)[2], length(tm.c), 1))
>>
>> for (i in 1:10) {
>> z[,,i,] <- i * ( sin(sqrt(x^2+y^2)))
>> }
>>
>> sin3dA <- new("Spatial3dArray",
>> ? ? ? data = z,
>> ? ? ? coords = list(x, y),
>> ? ? ? bbox = matrix(c(min(x), min(y), max(x), max(y), 2, 2), 2, 2,
>> ? ? ? dimnames = list(NULL, c("min","max"))),
>> ? ? ? time = tm.c,
>> ? ? ? btime = c(min(tm.c), max(tm.c)))
>>
>> dimnames(slot(sin3dA, "data")) = list(NULL,
>> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? NULL,
>> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? slot(sin3dA, "time"),
>> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? c("a"))
>> names(slot(sin3dA, "coords")) <- c("x", "y")
>>
>>
>> ##########################################
>>
>> for the coordinates method I would like to have two options on how to return
>> the coordinates; "list" or default "sp":
>>
>> coordinates.3dArray <- function (obj, type = "sp") {
>> ? ? ? lat <- slot(obj, "coords")[[1]]
>> ? ? ? long <- slot(obj, "coords")[[2]]
>> ? ? ? if (type == "list") {
>> ? ? ? ? ? ? ? return(list(long=long, lat=lat))
>> ? ? ? } else if (type == "sp") {
>> ? ? ? ? ? ? ? res <- as.matrix(cbind(c(long), c(lat)))
>> ? ? ? ? ? ? ? dimnames(res) <- list(NULL, c("x1", "x2"))
>> ? ? ? }
>> }
>> setMethod("coordinates", signature("Spatial3dArray"), coordinates.3dArray)
>>
>> This means that the default coordinates method in sp has to include the option
>> "..." . Would it be possible to include this in a future release of sp?
>>
>> The reason I want to keep the list option is to use a matrix oriented approach
>> in spplot, overlay, etc. methods. I also feel having a matrix/array approach
>> with these kind of data makes sense. Allowing type = "sp" means overlay() will
>> work more or less out of the box (however I would like to return a matrix),
>> and still I could get the list/matrix when desired.
>>
>>
>> Best wishes
>> Torleif Markussen Lunde
>> Centre for International Health
>> Bjerknes Centre for Climate Research
>> University of Bergen
>> Norway
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of M?nster
> Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
> http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
------------------------------
Message: 18
Date: Thu, 19 Nov 2009 23:34:37 +0100
From: Torleif Markussen Lunde <torleif.lunde at cih.uib.no>
Subject: Re: [R-sig-Geo] Spatial3dArray - coordinates method
To: Michael Sumner <mdsumner at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <200911192334.37645.torleif.lunde at cih.uib.no>
Content-Type: Text/Plain; charset="utf-8"
Hi Mike
At the moment I have written wrapper functions around ncdf. As in get.var.ncdf
you can subset which area to read, and which (continuous) time dimensions you
want to read. At the moment the functions (to correctly return time, lat and
long) are limited to output from the WRF-model (http://www.wrf-model.org/),
but it could easily be modified to other netcdf files as long you know the name
of lat, long and time.
If r-forge accepts my application (delivered yesterday) the methods should
appear there soon as rwrf.
Since ncdf works "out of the box" on Fedora I landed on that one. I also
tested it on windows XP, and no problems there either. It is a couple of years
since I tried RNetCDF, so I should not speak to strongly in favor for any of
them. Of course the bad thing with both is that they have not been updated for
a while. Still my impression is that they are both superior to GDALs NetCDF
support(?).
If the classes proves to be robust and the quality is sufficient it would be
natural to include them in the sp package in the future (instead of having to
dig around to look for sp classes). Currently spplot methods for
Spatial3dArray only support plotting of single times for levelplot,
contourplot, and wireframe (unless animation is requested). I have some
problems with the lattice graphics inside functions (also with print()) that
needs to be sorted out first.
Examples:
first convert Spatial3dArray to data.frame
print(levelplot(z‾long*lat | time, data = tmp.df, aspect="iso"))
will cause the last time to over plot all frames inside a function.
What could work is:
for (i in 1:2) {
if (i == 1) {
a <- paste('c(slot(surf, "data")[ , , ', i, ',2])',
sep = "")
} else {
a <- paste(a, ' + c(slot(surf, "data")[ , , ', i, ',2])',
sep = "")
}
}
levelplot(parse(text=a)‾
c(slot(surf, "coords")$long) * c(slot(surf, "coords")$lat),
strip = strip.custom(factor.levels=c("time1", "time2")))
where neither parse, print(a, quote = FALSE) or cat works as the z variable.
Probably this approach is not meant to work.
There are quite some issues that has to be solved before the class is
production-ready.
By the way, if someone has an idea on how to solve the last code snippets I
would be more than happy.
Best wishes
Torleif
On Thursday 19 November 2009 20:23:35 Michael Sumner wrote:
> Wow, this is great - I was thinking about this just yesterday.
>
> Torleif: do you have an opinion on which NetCDF path is the most
> useful for R with sp? RNetCDF or ncdf? GDAL is workable but takes
> extra effort to build and then reconstruct 3d/4d from 2d bands. (I use
> Windows mostly)
>
> I use the RNetCDF package a lot, mainly because that's the one I
> learnt to use first - there are binaries for Windows. It has some
> problems in terms of R-style but they could be easily fixed.
>
> Regards, Mike.
>
> On Fri, Nov 20, 2009 at 12:26 AM, Edzer Pebesma
>
> <edzer.pebesma at uni-muenster.de> wrote:
> > No problems; sp in csv now has this, the next release will have it.
> >
> > Torleif Markussen Lunde wrote:
> >> Hi
> >>
> >> To read netcdf data (or any other "gridded" spatial time data) I find it
> >> convenient to define new classes Spatial3dArray and Spatial4dArray.
> >>
> >> setClass("Spatial3dArray",
> >> representation("Spatial", data = "array", coords = "list",
> >> time = "character", btime = "character"),
> >> prototype= list(data = array(NA, c(1,1,1,1)),
> >> bbox=matrix(NA),
> >> proj4string = CRS(as.character(NA)),
> >> coords = list(1,1),
> >> time = "posix",
> >> btime = "posix"))
> >>
> >>
> >> ##########################################
> >> ###################EXAMPLE##################
> >> ##########################################
> >>
> >> x <- matrix(seq(-10, 10, length = 100), 100, 100,
> >> byrow = FALSE)
> >> y <- matrix(seq(-10, 10, length = 100), 100, 100,
> >> byrow = TRUE)
> >>
> >> tm <- 1:10
> >> tm.c <- as.character(seq(as.POSIXct("2002-01-01 06:00:00",
> >> "2002-01-01 15:00:00"),
> >> by="hours",
> >> length.out=10))
> >>
> >> z <- array(NA, c(dim(x)[1], dim(x)[2], length(tm.c), 1))
> >>
> >> for (i in 1:10) {
> >> z[,,i,] <- i * ( sin(sqrt(x^2+y^2)))
> >> }
> >>
> >> sin3dA <- new("Spatial3dArray",
> >> data = z,
> >> coords = list(x, y),
> >> bbox = matrix(c(min(x), min(y), max(x), max(y), 2, 2), 2, 2,
> >> dimnames = list(NULL, c("min","max"))),
> >> time = tm.c,
> >> btime = c(min(tm.c), max(tm.c)))
> >>
> >> dimnames(slot(sin3dA, "data")) = list(NULL,
> >> NULL,
> >> slot(sin3dA, "time"),
> >> c("a"))
> >> names(slot(sin3dA, "coords")) <- c("x", "y")
> >>
> >>
> >> ##########################################
> >>
> >> for the coordinates method I would like to have two options on how to
> >> return the coordinates; "list" or default "sp":
> >>
> >> coordinates.3dArray <- function (obj, type = "sp") {
> >> lat <- slot(obj, "coords")[[1]]
> >> long <- slot(obj, "coords")[[2]]
> >> if (type == "list") {
> >> return(list(long=long, lat=lat))
> >> } else if (type == "sp") {
> >> res <- as.matrix(cbind(c(long), c(lat)))
> >> dimnames(res) <- list(NULL, c("x1", "x2"))
> >> }
> >> }
> >> setMethod("coordinates", signature("Spatial3dArray"),
> >> coordinates.3dArray)
> >>
> >> This means that the default coordinates method in sp has to include the
> >> option "..." . Would it be possible to include this in a future release
> >> of sp?
> >>
> >> The reason I want to keep the list option is to use a matrix oriented
> >> approach in spplot, overlay, etc. methods. I also feel having a
> >> matrix/array approach with these kind of data makes sense. Allowing type
> >> = "sp" means overlay() will work more or less out of the box (however I
> >> would like to return a matrix), and still I could get the list/matrix
> >> when desired.
> >>
> >>
> >> Best wishes
> >> Torleif Markussen Lunde
> >> Centre for International Health
> >> Bjerknes Centre for Climate Research
> >> University of Bergen
> >> Norway
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> > --
> > Edzer Pebesma
> > Institute for Geoinformatics (ifgi), University of M?nster
> > Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
> > 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
> > http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
------------------------------
Message: 19
Date: Thu, 19 Nov 2009 22:18:42 -0500
From: Ned Horning <horning at amnh.org>
Subject: [R-sig-Geo] Problem using mahasuhab
To: r-sig-geo at stat.math.ethz.ch
Message-ID: <4B060A92.10103 at amnh.org>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi - I was wondering if anyone out there can help with my effort to
create habitat suitability maps using mahasuhab from the adehabitat
package or another package if there is a better option. I would like to
compare the resulting maps with some software a colleague is working on.
When I try to import an ascii grid using the import.asc method I get the
following error:
--
Error in if ((yll[[1]][1] == "yllcenter") | (xll[[1]][1] ==
"YLLCENTER")) corny <- FALSE :
missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In import.asc(filename, type = "numeric") : NAs introduced by coercion
2: In import.asc(filename, type = "numeric") : NAs introduced by coercion
--
I am able to read the test ascii grid file that comes with the package
just fine. I can also read my file using stack() from the raster package
but I don't know how to convert that RasterStack to an asc or kasc
object which seems to be necessary to run mahasuhab.
Any pointers?
Ned
------------------------------
Message: 20
Date: Thu, 19 Nov 2009 23:11:48 -0500
From: Mathieu Basille <basille at biomserv.univ-lyon1.fr>
Subject: Re: [R-sig-Geo] Problem using mahasuhab
To: Ned Horning <horning at amnh.org>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4B061704.10500 at biomserv.univ-lyon1.fr>
Content-Type: text/plain; charset=ISO-8859-1
Hi Ned,
Adehabitat was not meant to import all kind of raster maps. It might be
that your file is not in the right format. Are your file with such a header:
ncols 308
nrows 435
xllcorner 400000
yllcorner 6515000
cellsize 1000
NODATA_value -9999
(whatever the values are)
If yes, then it might be another deeper problem... Can you read it
through rgdal? If yes, you can use it and then convert it to a kasc
class, following this example:
bla <- readGDAL("your_map.asc")
blaK <- spixdf2kasc(bla)
Hope this helps,
Mathieu.
Ned Horning a ?crit :
> Hi - I was wondering if anyone out there can help with my effort to
> create habitat suitability maps using mahasuhab from the adehabitat
> package or another package if there is a better option. I would like to
> compare the resulting maps with some software a colleague is working on.
>
> When I try to import an ascii grid using the import.asc method I get the
> following error:
> --
> Error in if ((yll[[1]][1] == "yllcenter") | (xll[[1]][1] ==
> "YLLCENTER")) corny <- FALSE :
> missing value where TRUE/FALSE needed
> In addition: Warning messages:
> 1: In import.asc(filename, type = "numeric") : NAs introduced by coercion
> 2: In import.asc(filename, type = "numeric") : NAs introduced by coercion
> --
>
> I am able to read the test ascii grid file that comes with the package
> just fine. I can also read my file using stack() from the raster package
> but I don't know how to convert that RasterStack to an asc or kasc
> object which seems to be necessary to run mahasuhab.
>
> Any pointers?
>
> Ned
--
‾$ whoami
Mathieu Basille, Post-Doc
‾$ locate
Laboratoire d'?cologie Comportementale et de Conservation de la Faune
+ Centre d'?tude de la For?t
D?partement de Biologie
Universit? Laval, Qu?bec
‾$ info
http://ase-research.org/basille
‾$ fortune
``If you can't win by reason, go for volume.''
Calvin, by Bill Watterson.
------------------------------
Message: 21
Date: Thu, 19 Nov 2009 23:22:34 -0500
From: Ned Horning <horning at amnh.org>
Subject: Re: [R-sig-Geo] Problem using mahasuhab
To: Mathieu Basille <basille at biomserv.univ-lyon1.fr>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4B06198A.5040903 at amnh.org>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Mathieu,
Thanks for the help - rgdal / spixdf2kasc seems to have done the trick.
Here is the header from my ascii grid:
--
ncols 900
nrows 1650
xllcorner 0
yllcorner 0
cellsize 1000
NODATA_value -9999
--
All the best,
Ned
Mathieu Basille wrote:
> Hi Ned,
>
> Adehabitat was not meant to import all kind of raster maps. It might be
> that your file is not in the right format. Are your file with such a header:
>
> ncols 308
> nrows 435
> xllcorner 400000
> yllcorner 6515000
> cellsize 1000
> NODATA_value -9999
>
> (whatever the values are)
>
> If yes, then it might be another deeper problem... Can you read it
> through rgdal? If yes, you can use it and then convert it to a kasc
> class, following this example:
>
> bla <- readGDAL("your_map.asc")
> blaK <- spixdf2kasc(bla)
>
> Hope this helps,
> Mathieu.
>
>
> Ned Horning a ?crit :
>
>> Hi - I was wondering if anyone out there can help with my effort to
>> create habitat suitability maps using mahasuhab from the adehabitat
>> package or another package if there is a better option. I would like to
>> compare the resulting maps with some software a colleague is working on.
>>
>> When I try to import an ascii grid using the import.asc method I get the
>> following error:
>> --
>> Error in if ((yll[[1]][1] == "yllcenter") | (xll[[1]][1] ==
>> "YLLCENTER")) corny <- FALSE :
>> missing value where TRUE/FALSE needed
>> In addition: Warning messages:
>> 1: In import.asc(filename, type = "numeric") : NAs introduced by coercion
>> 2: In import.asc(filename, type = "numeric") : NAs introduced by coercion
>> --
>>
>> I am able to read the test ascii grid file that comes with the package
>> just fine. I can also read my file using stack() from the raster package
>> but I don't know how to convert that RasterStack to an asc or kasc
>> object which seems to be necessary to run mahasuhab.
>>
>> Any pointers?
>>
>> Ned
>>
>
>
>
------------------------------
Message: 22
Date: Thu, 19 Nov 2009 21:54:42 -0700
From: Alan Swanson <alan.swanson at umontana.edu>
Subject: [R-sig-Geo] changing the data type of a gdal dataset
To: r-sig-geo at stat.math.ethz.ch
Message-ID: <4B062112.3080809 at umontana.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Dear R gurus,
I have a function that applies various model prediction functions over a
set of large image files, producing a single output file with the same
spatial extent. Due to memory issues, I'm breaking the input and output
files into tiles. I have this working except for one small issue
regarding data types.
I create a new gdal transient dataset by copying an existing one using:
handle <- GDAL.open(fullnames[1],read.only=T)
tds <-
copyDataset(handle,driver=new('GDALDriver','GTiff'),strict=F,options=NULL)
...
putRasterData(tds,t(preds), offset= c(strt[1], 0))
...
saveDataset(tds,outfile.p)
GDAL.close(tds)
Which works great, except that my output always needs to be floating
point, but the input may be byte or integer, in which case the output
dataset retains the format of the input file. So I either need to
change the data type of the new file, or create the new file using:
tds <- new("GDALTransientDataset",driver,dims[1],dims[2], type="Float32")
and then copy the spatial reference information from an existing
dataset. I can't figure out how to do either of these. Your help would
be much appreciated.
Cheers,
Alan
------------------------------
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
End of R-sig-Geo Digest, Vol 75, Issue 18
More information about the R-sig-Geo
mailing list