[R-sig-Geo] R-stats - GRASS6
Markus Neteler
neteler at itc.it
Thu Jun 9 17:14:54 CEST 2005
Dear Roger,
sorry for the delay (was out for teaching). My installation is
untouched since then.
I tested you functions as proposed below, it helped:
> library(spgrass6)
Loading required package: rgdal
Loading required package: abind
Loading required package: pixmap
Geospatial Data Abstraction Library extensions to R successfully loaded
Loading required package: sp
> G <- gmeta6()
[now paste of the modified readCELL6sp() function]
> landcover <- readCELL6sp("landcover.30m", cat = TRUE)
> summary(landcover)
Object of class SpatialGridDataFrame
Coordinates:
min max
coords.x1 589995 608985
coords.x2 4913715 4927995
Is projected: TRUE
proj4string : [+proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +no_defs +nadgrids=conus]
Number of points: 2
Grid attributes:
cellcentre.offset cellsize cells.dim
1 589995 30 634
2 4913715 30 477
Data attributes:
landcover.30m
Evergreen Forest :135375
Grasslands/Herbaceous : 67396
Pasture/Hay : 56263
Small Grains : 9023
Emergent Herbaceous Wetlands: 5673
(Other) : 18587
NA's : 10101
Now I tried to plot the map according your new article in the
GRASS Newsletter (http://grass.itc.it/newsletter/):
> image(landcover,"landcover.30m")
Error in image.default(as.image.SpatialGridDataFrame(x[attr], xcol, ycol), :
invalid z limits
In addition: Warning messages:
1: no finite arguments to min; returning Inf
2: no finite arguments to max; returning -Inf
Maybe I am missing something?
Thanks
Markus
On Wed, Jun 01, 2005 at 01:33:26PM +0200, Roger Bivand wrote:
> On Tue, 31 May 2005, Markus Neteler wrote:
>
> > Dear Roger,
> >
> > On Tue, May 31, 2005 at 06:53:31PM +0200, Roger Bivand wrote:
> > > On Tue, 31 May 2005, Markus Neteler wrote:
> > >
> > > > Dear list members,
> > > >
> > > > I am currently trying out the new GRASS 6/R-stats interface.
> > > > I installed all software today (R-2.1.0, all interfaces, using
> > > > GRASS 6.1-CVS).
> > > >
> > > > There is a little problem which I don't understand:
> > > >
> > > > > library(spgrass6)
> > > > > G <- gmeta6()
> > > > > str(G)
> > > > List of 28
> > > > $ GISDBASE : chr "/ssi0/ssi/neteler/grassdata"
> > > > $ LOCATION_NAME : chr "spearfish60"
> > > > $ MAPSET : chr "neteler"
> > > > $ GRASS_DB_ENCODING: chr "utf-8"
> > > > $ DEBUG : chr "0"
> > > > $ MONITOR : chr "x0"
> > > > $ GRASS_GUI : chr "text"
> > > > $ projection : chr "1 (UTM)"
> > > > [...]
> > > >
> > > > > system("r.info -t landcover.orig")
> > > > datatype=CELL
> > > > > landcover <- readCELL6sp("landcover.orig", cat = TRUE)
> > > > Error in "names<-.default"(`*tmp*`, value = "landcover.orig") :
> > > > 'names' attribute [1] must be the same length as the vector [0]
> > >
> > > I can't replicate this, I'm afraid (R 2.1.0, Linux, sp 0.7-7, spgrass6
> > > 0.1-4 (from CVS on sourceforge, GRASS 6.0.0 release). I don't think the sp
> > > or spgrass6 versions make a difference.
> > >
> > > Could you say debug(readCELL6sp) and re-run this please - that should show
> > > which line is failing. Once we know that, please check the arguments to
> > > that function for sanity. I don't think it is important, but I'm running
> > > under LANG=en_GB ; export LANG, not UTF8. My current guess is that the
> > > problem is in:
> > >
> > > res <- read.asciigrid(tmpfl, colname = vname, proj4string = p4)
> > >
> > > where vname is the layer name. If you additionally do
> > > debug(read.asciigrid), it may fail at: names(df) = colname - please see
> > > what is in colname, and if length(colname) is the same as length(df). One
> > > possibility is that there is something amiss in read.asciigrid().
> > >
> > > The landcover.orig layer is strange because there are no category labels
> > > at all, so:
> >
> > [... took the wrong landuse map, but it doesn't change ...]
> >
> > > library(spgrass6)
> > Loading required package: rgdal
> > Loading required package: abind
> > Loading required package: pixmap
> > Geospatial Data Abstraction Library extensions to R successfully loaded
> > Loading required package: sp
> > > G <- gmeta6()
> >
> > > landcover <- readCELL6sp("landcover.30m", cat = TRUE)
> > Error in "names<-.default"(`*tmp*`, value = "landcover.30m") :
> > 'names' attribute [1] must be the same length as the vector [0]
> > > debug(readCELL6sp)
> >
> > > system("r.stats -l -q landcover.30m | head -20")
> > 11 Open Water
> > 21 Low Intensity Residential
> > 22 High Intensity Residential
> > 23 Commercial/Industrial/Transportation
> > 31 Bare Rock/Sand/Clay
> > 32 Quarries/Strip Mines/Gravel Pits
> > 41 Deciduous Forest
> > 42 Evergreen Forest
> > 43 Mixed Forest
> > 51 Shrubland
> > 71 Grasslands/Herbaceous
> > 81 Pasture/Hay
> > 82 Row Crops
> > 83 Small Grains
> > 85 Urban/Recreational Grasses
> > 91 Woody Wetlands
> > 92 Emergent Herbaceous Wetlands
> > * no data
> >
> > > landcover <- readCELL6sp("landcover.30m", cat = TRUE)
> > debugging in: readCELL6sp("landcover.30m", cat = TRUE)
> > debug: {
> > [...]
> > }
> > Browse[1]>
> > debug: tmpfl <- tempfile()
> > Browse[1]>
> > debug: system(paste("r.out.arc input=", vname, " output=", tmpfl, sep = ""))
> > Browse[1]>
> > debug: library(sp)
> > Browse[1]>
> > debug: p4 <- CRS(system("g.proj -j -f", intern = TRUE))
> > Browse[1]>
> > debug: res <- read.asciigrid(tmpfl, colname = vname, proj4string = p4)
> > Browse[1]>
> > debug: if (cat) {
> > cats <- strsplit(system(paste("r.stats -l -q", vname), intern = TRUE),
> > " ")
> > catnos <- sapply(cats, function(x) x[1])
> > catlabs <- sapply(cats, function(x) paste(x[-1], collapse = " "))
> > if (any(!is.na(match(catnos, "*")))) {
> > isNA <- which(catnos == "*")
> > catnos <- catnos[-isNA]
> > catlabs <- catlabs[-isNA]
> > }
> > res at data[, 1] <- factor(res at data[, 1], levels = catnos, labels = catlabs)
> > } else {
> > res at data[, 1] <- as.integer(res at data[, 1])
> > }
> > Browse[1]>
> > debug: cats <- strsplit(system(paste("r.stats -l -q", vname), intern = TRUE),
> > " ")
> > Browse[1]>
> > debug: catnos <- sapply(cats, function(x) x[1])
> > Browse[1]>
> > debug: catlabs <- sapply(cats, function(x) paste(x[-1], collapse = " "))
> > Browse[1]>
> > debug: if (any(!is.na(match(catnos, "*")))) {
> > isNA <- which(catnos == "*")
> > catnos <- catnos[-isNA]
> > catlabs <- catlabs[-isNA]
> > }
> > Browse[1]>
> > debug: isNA <- which(catnos == "*")
> > Browse[1]>
> > debug: catnos <- catnos[-isNA]
> > Browse[1]>
> > debug: catlabs <- catlabs[-isNA]
> > Browse[1]>
> > debug: res at data[, 1] <- factor(res at data[, 1], levels = catnos, labels = catlabs)
> > Browse[1]>
> > Error in "names<-.default"(`*tmp*`, value = "landcover.30m") :
> > 'names' attribute [1] must be the same length as the vector [0]
> > >
> >
> >
> > > debug(read.asciigrid)
> > > undebug(readCELL6sp)
> >
> > > landcover <- readCELL6sp("landcover.30m", cat = TRUE)
> > [...]
> > Browse[1]>
> > debug: df = data.frame(map)
> > Browse[1]>
> > debug: names(df) = colname
> > Browse[1]> length(colname)
> > [1] 1
> > Browse[1]> length(df)
> > [1] 1
> > Browse[1]> colname
> > [1] "landcover.30m"
> > Browse[1]> df
> > [...all values...]
> >
> > > str(df)
> > `data.frame': 302418 obs. of 1 variable:
> > $ map: num NA NA NA NA NA NA NA NA NA NA ...
> > Browse[1]>
> > debug: grid = GridTopology(c(xllcenter, yllcenter), rep(cellsize, 2),
> > c(ncols, nrows))
> > Browse[1]>
> > debug: SpatialGridDataFrame(grid, data = df, proj4string = proj4string)
> > Browse[1]>
> > exiting from: read.asciigrid(tmpfl, colname = vname, proj4string = p4)
> > Error in "names<-.default"(`*tmp*`, value = "landcover.30m") :
> > 'names' attribute [1] must be the same length as the vector [0]
> > >
> >
> > I don't know if this is already helpful.
>
> Yes, looks as though the changes we made to sp in Valencia (you have sp
> 0.7-7, right?) are getting in the way. When the data is being made into a
> factor, the name of the data seems to be getting messed up. I have:
>
> readCELL6sp <- function(vname, cat=FALSE) {
> tmpfl <- tempfile()
> system(paste("r.out.arc input=", vname, " output=", tmpfl, sep=""))
> library(sp)
> p4 <- CRS(system("g.proj -j -f", intern=TRUE))
> res <- read.asciigrid(tmpfl, colname=vname, proj4string=p4)
> if (cat) {
> cats <- strsplit(system(paste("r.stats -l -q", vname),
> intern=TRUE), " ")
> catnos <- sapply(cats, function(x) x[1])
> catlabs <- sapply(cats, function(x) paste(x[-1], collapse=" "))
> if (any(!is.na(match(catnos, "*")))) {
> isNA <- which(catnos == "*")
> catnos <- catnos[-isNA]
> catlabs <- catlabs[-isNA]
> }
> res at data[[1]] <- factor(res at data[[1]], levels=catnos, labels=catlabs)
> } else {
> res at data[[1]] <- as.integer(res at data[[1]])
> }
> res
> }
>
> on my machine, I think it is on CVS, but cannot check now. The change
> Edzer and I made to sp was the same that was made to the R/GRASS interface
> years ago, that is do not use data frames where possible, because they
> generate unique character row names, which we do not need. So we changed
> from the @data slot being a data frame that can be accessed like a matrix,
> to a list of equal length vectors, which can only be accessed by list
> member. Please try this change in the function. The error message isn't
> the same in the following, but you can see what is going on (unless I'm
> still wrong):
>
> > df <- data.frame(a=sample(1:5,100,replace=TRUE))
> > str(df)
> `data.frame': 100 obs. of 1 variable:
> $ a: int 3 2 2 2 3 1 3 1 5 2 ...
> > ldf <- as.list(df)
> > str(ldf)
> List of 1
> $ a: int [1:100] 3 2 2 2 3 1 3 1 5 2 ...
> > df1 <- df
> > df1[,1] <- factor(df[,1])
> > str(df1)
> `data.frame': 100 obs. of 1 variable:
> $ a: Factor w/ 5 levels "1","2","3","4",..: 3 2 2 2 3 1 3 1 5 2 ...
> > ldf1 <- ldf
> > ldf1[,1] <- factor(ldf[,1])
> Error in ldf[, 1] : incorrect number of dimensions
> > ldf1[[1]] <- factor(ldf[[1]])
> > str(ldf1)
> List of 1
> $ a: Factor w/ 5 levels "1","2","3","4",..: 3 2 2 2 3 1 3 1 5 2 ...
>
> Best wishes,
>
> Roger
>
> PS. I think 0.7-7 read.asciigrid() is also wrong, that is it should not
> use data.frame() anyway, but I can't check that now.
>
> >
> > Best regards
> >
> > Markus
> >
>
> --
> 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
>
--
Markus Neteler <neteler itc it> http://mpa.itc.it
ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
MPBA - Predictive Models for Biol. & Environ. Data Analysis
Via Sommarive, 18 - 38050 Povo (Trento), Italy
More information about the R-sig-Geo
mailing list