[R] can't vectorize an expression

[Roger Bivand]

On Wed, 12 Dec 2001, Robin Hankin wrote:

> Dear R support network
> 
> I have a problem that is driving me crazy.
> 
> I have a dataframe with about 74000 landscapes which I call "land".  A
> landscape is a 2km -by- 2km square.
> 
> Land has three columns: land$lat, land$long, and land$description.
> The last one holds a NON-unique (integer) description of each
> landscape.  There are maybe 100 distinct descriptions.  Identifying
> landscapes that have identical descriptions gives equivalence classes.
> 
> I am interested in the largest few equivalence classes (that is, a
> *large* set of landscapes with identical descriptions).  The largest
> equivalence class holds around 1500 landscapes.
> 
> 
> Thus:
> 
> top <- rev(sort(table(land$description)))
> top.values  <- as.numeric(names(top))
> 
> So far so good: top.values holds the descriptions of the largest
> equivalence classes first.
> 

Using a similar setup, for land cover classes for a raster map layer:

> table(z)
z
    1     2     3     4     5     6     7     8 
 1671  6632   555  1749 22913 16816  6558   706 
> top <- rev(sort(table(z)))
> top
    5     6     2     7     4     1     8     3 
22913 16816  6632  6558  1749  1671   706   555 
> top.values  <- as.integer(names(top))
> ztop <- apply(matrix(z, ncol=1), 1, function(x) which(top.values == x))

Apply takes a bit of time, but certainly less than keying in all the
commands. To drop the classes over the sixth:

> ztop[ztop > 6] <- NA
> summary(as.ordered(ztop))
    1     2     3     4     5     6  NA's 
22913 16816  6632  6558  1749  1671  1261 

Since I'm using the GRASS GIS interface here (in Devel), I can retrieve
the position data:

> summary(G)
Data from GRASS 5.0 LOCATION leics with 240 columns and 240 rows;
The west-east range is: 444000, 456000, and the south-north: 310000,
322000;
West-east cell sizes are 50 units, and south-north 50 units.

Given that, and knowing the sequences of cell centres in x and y:

> image(G$xseq, G$yseq, t(matrix(ztop, nrow=240, ncol=240,
byrow=T))[,240:1], asp=1, breaks=seq(0.5,6.5,1),
col=c("red","blue","yellow","green","purple","cyan")

gives the plot (asp=1 to keep aspect ratio). Looking at the code in the
pixmap package, you'll also see some ways to get the image out.


> Now things turn to custard; the following lines show the only way I
> could think of to extract the six largest equivalence classes:
> 
> which(land$description == top.values[1]) -> biggest1
> which(land$description == top.values[2]) -> biggest2
> which(land$description == top.values[3]) -> biggest3
> which(land$description == top.values[4]) -> biggest4
> which(land$description == top.values[5]) -> biggest5
> which(land$description == top.values[6]) -> biggest6 
> 
> plotting them is relatively easy:
> 
> plot(b[,1:2],pch=16,cex=0.2)
> matplot(land$lat[biggest1],land$long[biggest1], col="red",   type="p",pch=16,add=TRUE)
> matplot(land$lat[biggest2],land$long[biggest2], col="blue",  type="p",pch=16,add=TRUE)
> matplot(land$lat[biggest3],land$long[biggest3], col="yellow",type="p",pch=16,add=TRUE)
> matplot(land$lat[biggest4],land$long[biggest4], col="green", type="p",pch=16,add=TRUE)
> matplot(land$lat[biggest5],land$long[biggest5], col="purple",type="p",pch=16,add=TRUE)
> matplot(land$lat[biggest6],land$long[biggest6], col="cyan",  type="p",pch=16,add=TRUE)
> 
> 
> QUESTION: How do I vectorize this?
> 
> --
> 
> 
> A test dataset might look like
> 
> 
> lat long description
>  172.1100 -34.1500 8
>  172.1300 -34.1500 8
>  172.1500 -34.1500 8
>  172.1300 -34.1700 8
>  172.1500 -34.1700 7
>  173.0100 -34.3900 7
>  173.0300 -34.3900 6
>  172.8700 -34.4100 7
>  172.8900 -34.4100 7
>  172.9700 -34.4100 6
>  172.9900 -34.4100 8
>  173.0100 -34.4100 8
>  173.0300 -34.4100 8
>  173.0500 -34.4100 8
>  172.6700 -34.4300 8
>  172.6900 -34.4300 6
>  172.7100 -34.4300 6
>  172.7300 -34.4300 6
>  172.7500 -34.4300 6
> 		   
> 
> 

Roger

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no
and: Department of Geography and Regional Development, University of
Gdansk, al. Mar. J. Pilsudskiego 46, PL-81 378 Gdynia, Poland.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._