Dear Zhijie Zhang, I'm working on an approach to this. I was asked (also on Tuesday) whether it is possible to disguise map position to increase confidentiality for patient locations. This is a valid reason to elide spatial position, so I'll try to see what can be done, and it will include your case.
I will send a draft version for you to try when I have something that seems to work. Best wishes, Roger On Tue, 25 Sep 2007, zhijie zhang wrote: > Dear Porf. Roger, > There are still a little problem with function mess_up(). > #Programs to read the rivers > library(sp);library(foreign);library(maptools) > > rivers<- readShapeLines("e:/qiupu.shp") > > #loops > > lns <- slot(rivers, "lines") > > new_lns <- lapply(lns, function(x) > { > Lns <- slot(x, "Lines") > new_Lns <- lapply(Lns, function(y) > { > new_crds <-* mess_up*(slot(y, "coords")) > Line(new_crds) > }) > Lines(new_Lns, ID=slot(x, "ID")) > }) > new_rivers <- SpatialLines(new_lns) > > 1. #for your mess_up() - note, untried. > *#Erros: no function mess_up()* > > 2. #Also note that the projection will not be defined. > Actually, i have projected in ARCGIS, so the projection will be unneccessary > in R. > > Thanks. > > > > > > On 9/24/07, Roger Bivand <[EMAIL PROTECTED]> wrote: >> >> On Mon, 24 Sep 2007, zhijie zhang wrote: >> >>> Dear Roger, >>> The reason for doing this is to get more general programs, which can be >>> used in the other similar conditions. If we can convert all the related >>> objects into internal [0,1], we can more easily use the intermediate >> results >>> calculated from unit square,such as .bandwidth in kernels. >> >> Exactly - with the original data, the units *are* general, because they >> are in meters. Provided the other data is also using the same metric, this >> also permits a direct interpretation of bandwidth, possibly with a >> physical or biological basis. In the same sense, you can also get a >> relative bandwidth from (bw / max(diff(bbox(<obj>)))), and multiply it out >> again for a different context. >> >>> I think the main problem maybe the lines object, that is qiupu.shp. The >>> scaling factor on the longer axis is from the guichi polygon. >>> #read the lines of two rivers in the guichi >>> rivers <- readShapeLines("e:/qiupu.shp") #change the >> coordinates??difficult >>> #plot(rivers,add=T) >>> I'm not very sure if they,especially qiupu lines, can be successfullyy >>> done. >> >> Just loops of lapply(), not worse: >> >> lns <- slot(rivers, "lines") >> >> new_lns <- lapply(lns, function(x) { >> Lns <- slot(x, "Lines") >> new_Lns <- lapply(Lns, function(y) { >> new_crds <- mess_up(slot(y, "coords")) >> Line(new_crds) >> }) >> Lines(new_Lns, ID=slot(x, "ID") >> }) >> new_rivers <- SpatialLines(new_lns) >> >> for your mess_up() - note, untried. Also note that the projection will not >> be defined. >> >> Roger >> >>> Thanks very much. >>> >>> On 9/24/07, Roger Bivand <[EMAIL PROTECTED]> wrote: >>>> >>>> Dear Zhijie Zhang, >>>> >>>> On Sun, 23 Sep 2007, zhijie zhang wrote: >>>> >>>>> Dear Roger, >>>>> I have a problem on converting the coordinates into interval [0,1], >>>> hoping >>>>> that you can give me a hand. >>>> >>>> I think that you need to explain why. If it is just because the >>>> coordinates "explode" because they are large, then simply scaling like >>>> this: >>>> >>>> library(rgdal) >>>> guichi <- readOGR(".", "guichi") >>>> plot(guichi, axes=TRUE) >>>> bbox(guichi) >>>> proj4string(guichi) >>>> t1 <- paste("+proj=tmerc +lat_0=0 +lon_0=117 +k=1.000000 +x_0=0", >>>> "+y_0=-3348000 +a=6378140 +b=6356755.288157528 +units=km") >>>> t2 <- spTransform(guichi, CRS(t1)) >>>> plot(t2, axes=TRUE) >>>> bbox(t2) >>>> # min max >>>> # r1 10.148541 79.88203 >>>> # r2 0.834946 61.90978 >>>> >>>> gets down to 70 by 80, which is not such a numerical risk (if need be, >> use >>>> +units=kmi to get nautical miles, which reduces the numbers even more). >> It >>>> has the key advantage that you can always get back to the true >> coordinates >>>> by transforming back. >>>> >>>> In the internal code, you will see that many packages guard by using >> [-1, >>>> +1] on both dimensions. In your case and if you know that you need >> this, >>>> you would have to work out the scaling factor on the longer axis, so >> that >>>> points on the shorter axis would be scaled correctly. But if you want >> to >>>> keep the dimensions' aspect, you could just shift the origin and the >> units >>>> as I show above. >>>> >>>> Hope this helps, >>>> >>>> Roger >>>> >>>>> My data consist of points(cases), lines(rivers) and polygon (guichi >>>> city), >>>>> and i want to change their coordinates into interval [0,1]. I have >> put >>>> the >>>>> data in the attachment, so that you can use it. >>>>> The following is the programs to read the data and one possible method >>>> for >>>>> conversion: >>>>> >>>>> library(sp) >>>>> library(foreign) >>>>> library(mgcv) >>>>> library(maptools) >>>>> >>>>> #read the polygon containing the studied points >>>>> guichi <- readShapePoly("e:/guichi.shp") #boundary polygons >> containing >>>> the >>>>> points >>>>> point_poly <- >>>>> >>>> >> getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(guichi)[[1]])[[1]]) >>>>> #get the coordinates of guichi >>>>> #plot(point_poly,xlab="x", ylab="y",type="l") >>>>> >>>>> #read the lines of two rivers in the guichi >>>>> rivers <- readShapeLines("e:/qiupu.shp") #change the >>>> coordinates??difficult >>>>> #plot(rivers,add=T) >>>>> >>>>> #read the points of cases and controls >>>>> case_control <- read.csv("e:/casecontrol.csv",sep=",", header=TRUE) >>>>> #plot(case_control$x,case_control$y) >>>>> >>>>> *#Plot the whole figure* >>>>> plot(point_poly,xlab="x", ylab="y",type="l") >>>>> plot(rivers,add=T) >>>>> points(case_control$x,case_control$y,add=T) >>>>> >>>>> >>>> >> ################################################################################### >>>>> #one of the possible methods to convert the x/y coordinates into >>>> interval >>>>> [0,1] >>>>> *# But it seems there is a problem: it convert the x/y coordinates >> into >>>>> interval [0,1], respectively. >>>>> # In my opinion, they should be expand or shrink according to the >> same >>>>> minimum/maximum value*. >>>>> st <- function(x)(x-min(x))/(max(x)-min(x)) >>>>> case_control[,c(8,9)] <- data.frame(lapply(case_control[,c(6,7)],st)) >>>>> >>>> >> ################################################################################### >>>>> *Q1. The main problem is on the rivers, which has multiple lines. >> There >>>> are >>>>> difficulties for conversion. >>>>> Q2. when i convert the x/y coordinates, i'm not very sure whether i >>>> should >>>>> use the same scale of conversion or not. * >>>>> Could you please help us to solve it? >>>>> Any suggestions/help are greatly appreciated. >>>>> >>>>> >>>>> >>>> >>>> -- >>>> 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: [EMAIL PROTECTED] >>>> >>> >>> >>> >>> >> >> -- >> 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: [EMAIL PROTECTED] >> > > > > -- 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: [EMAIL PROTECTED] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo