[R-sig-Geo] R-stats - GRASS6
Roger Bivand
Roger.Bivand at nhh.no
Thu Jun 9 23:32:56 CEST 2005
On Thu, 9 Jun 2005, Markus Neteler wrote:
> 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?
Yes, you are plotting a factor, which image() is seeing as a character
matrix, and so rightly protests when it tries to find its range. This is
OK:
> landcover1 <- readCELL6sp("landcover.30m", cat = FALSE)
> image(landcover1,"landcover.30m",
+ col=rainbow(nlevels(landcover at data[[1]])))
but a bit meaningless. I'm not sure how to recover the category label
matches to the GRASS color palette at the shell level - can r.colors
report the current match? Thanks for flagging this need for improvement.
Best wishes,
Roger
>
> 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
> >
>
>
--
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