[R-sig-Geo] azymuth error, polygons and Lisbon projection

[Roger Bivand]

On Thu, 26 Jul 2007, Marta Rufino wrote:

> Hello,
>
>
> I have four doubts:
>
> 1)
> why gzAzymuth is giving me Northern direction in this case, for example?

gzAzymuth() in the maptools package is for geographical coordinates, as 
its help page says, and yours seem to be projected.

>
> require(sp); require(rgdal)
>
> kk=data.frame(x=c(575552.8,576295.9), y=c(4103900,4103611))
> coordinates(kk)=~x+y
> plot(kk, typ="o", col=4)
> gzAzimuth(coordinates(kk)[1,], coordinates(kk)[2,])
>
> 2) how to cut a polygon (for example a border file) into chuncks? (i.e.
> separate areas, according to a line)

With difficulty, probably using functions in the gpclib package to 
intersect the polygon with a rectangle. There is code inside Rgshhs() in 
maptools that uses this approach - the difficulty is that the polygon 
border may cross the line several times, creating multiple included 
fragments.

>
> 3) how to create multiple polygons, from a grid file, according to a
> criteria. For example, :
>
> require(gstat)
> data(meuse)
> coordinates(meuse) = ~x+y
> data(meuse.grid)
> gridded(meuse.grid) = ~x+y
> m <- vgm(.59, "Sph", 874, .04) # ordinary kriging:
> x <- krige(log(zinc)~1, meuse, meuse.grid, model = m)
> spplot(x["var1.pred"], main = "ordinary kriging predictions")
> # how to create a polygon with the areas greater than 5...

t1 <- as.SpatialPolygons.SpatialPixels(x)
df <- slot(x, "data")
rownames(df) <- sapply(slot(t1, "polygons"), function(x) slot(x, "ID"))
t2 <- SpatialPolygonsDataFrame(t1, data=df)
library(maptools)
library(gpclib)
t3 <- unionSpatialPolygons(t2, t2$var1.pred > 5)
spl1 <- list("sp.polygons", t3, lwd=2, col="blue")
spplot(t2, "var1.pred", sp.layout=list(spl1))

>
>
> 4) Is it defined in proj4, the datum lisboa, Sistema of Hayford- Gauss
> Militar? and UTM, using zone 29N (29 is down here, for example...)?
> CRS("+proj=utm +zone=29 +datum=WGS84")

Maybe you should cross-check on Grids & Datums:

http://www.asprs.org/resources/GRIDS/

April 2002. It seems to have most of what you need to construct a proj4 
encoding for the Lisboa datum, Hayford-Gauss system maps, but you'll need 
to be careful, and check with other users. It looks very like 
"+init=epsg:20790" without a datum, but the G&D and EPSG reverse the false 
eastings and northings.

>
> How do we specify it?
>
>

Hope this helps,

Roger

> Thank you very much in advance,... any help will be much appreciated...
> Best regards,
> Marta
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

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