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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._