2014-02-15 49 views
-2

我有壓縮的形狀文件。 我想將xlables和ylables改爲等價的經度和緯度值,而不是Eastings和Northings。請建議我如何顯示經度和緯度來代替這個北部和東部。當我改變shapefile到WGS1984的投影,但是扭曲了這個區域,我不想這樣做。 我的代碼是在ggplot2中顯示Long Lat代替Easting和Northing

library(raster) 
library(RColorBrewer) 
library(ggplot2) 
library(ncdf) 
poly<-shapefile("Test_poly.shp", stringsAsFactors=FALSE, verbose=FALSE) 
ggplot(poly) + 
    aes(long,lat,group=group, fill = id) + 
    geom_polygon() + 
    geom_path(color="white") + 
    coord_equal() + 
    scale_fill_brewer("Bydeler") 
+1

你的鏈接是'Test_poly.shp',但shapefile是一組文件*,所以'Test_poly.dbf','Test_poly.prg','test_poly.shp'等等。你需要提供*所有這些*或你的代碼不會運行。 – jlhoward

+0

另外,查看'scale_x_continuous'來更改座標軸上的標籤。 –

回答

1

(感謝更新shape文件。不過,我也發現,我需要在我的系統和安裝,也許值得一提的rgdal [R包上凸出?或許這只是反映出一個事實,我不在R中做很多空間工作)

最終,我不認爲ggplot2有一個內置的系統,用於在軸中使用不同的映射單元。我在這裏介紹的是變通方法,有點難看,但應該起作用。

要獲得緯度和長軸,一個選項是將shapefile轉換爲lat和long直接。例如:

poly <- shapefile("Test_poly.shp", stringsAsFactors=FALSE, verbose=FALSE) 
poly <- spTransform(poly, CRS("+proj=longlat +ellps=GRS80")) 

ggplot(poly2) + 
    aes(long,lat,group=group, fill = id) + 
    ylim(c())+ 
    xlim(c(1324794,1341700)) + 
    geom_polygon() + 
    geom_path(color="white") + 
    coord_equal() + 
    scale_fill_brewer("Bydeler") 

正如您指出的那樣,這會扭曲地圖。我認爲這是因爲一個緯度單位與一個單位緯度之間的差異不同,但是你已經告訴ggplot使用相同的座標。

我能看到的最好的解決方法是弄清楚在這種情況下經緯度與緯度之比是多少,並用coord_fixed(ratio = X)(其中X是比率)而不是使用coord_equal()來指定。

您應該能夠通過打印poly對象和未轉化和轉化形狀的物體之間比較「程度」弄清楚這個比例:

> poly 
    class  : SpatialPolygonsDataFrame 
    features : 7 
    extent  : 1211747, 1466285, 728771.9, 999422.3 (xmin, xmax, ymin, ymax) 
    coord. ref. : +proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
    variables : 1 
    names  : Id 
    min values : 0 
    max values : 0 

另外,您可以使用原來的shape文件,然後使用手冊以獲得緯度和經度,例如,scale_y_continuous(breaks = c(700000, 800000, 900000,1000000), labels = c(51,52,53,54))

這些不是一個優雅的解決方案,但應該適用於製作演示文稿或出版物的目的。