[R-sig-Geo] problem on converting the coordinates into interval [0, 1]?

[Roger Bivand]

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 <Roger.Bivand at nhh.no> 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 <Roger.Bivand at nhh.no> 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: 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
>>
>
>
>
>

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