Как я могу добавить сетку широты и долготы на проецируемую карту?

РЕДАКТИРОВАТЬ: вопрос действительно относится к тому, как можно добавить линии сетки широты / долготы на проецируемую карту. Я изменил название на соответствие.

У меня есть несколько слоев в географических координатах. Я хочу построить их в проекции LCC, но иметь географическую (широту / долготу) сетку. Из mapproj я могу использовать map.grid() для добавления сетки с ограничениями, установленными аргументом lim. Требуется вектор или объект диапазона:

вектор из 4 чисел, определяющих пределы: c (долгое низкое, долгое высокое, латинское низкое, латинское высокое). lim также может быть списком с компонентом с именем range, например результатом map, из которого берутся пределы.

Я создаю свою карту, обрезая большой векторный слой с помощью многоугольника отсечения:

myPoly <- readOGR(dsn=".", layer="myPolygon")  # just a shapefile in geographic coords
library(raster)  # To convert an 'extent' object to a SpatialPolygons object
cp <- as(extent(146, 149, -39, -37.5), "SpatialPolygons")
proj4string(cp) <- CRS(proj4string(myPoly))  # copy from shapefile

# Transform and plot:
lcc <- CRS("+init=epsg:3111")
myPoly.proj <- spTransform(myPoly, lcc)
cp.proj <- spTransform(cp, lcc)  # transform the clip box
myPoly.proj.clip <- gIntersection(myPoly.proj, cp.proj, byid=TRUE)
plot(myPoly.proj.clip) 

# Then finally, add a lat/long grid:
map.grid(lim=as.vector(cp.proj@bbox), labels=TRUE)

Эта последняя строка неверна, поскольку возвращаемый @bbox - это xmin, ymin, xmax, ymax, но должен быть в xmin, xmax, ymin, ymax. У всего этого должно быть простое решение, но, как обычно, я потерялся в водовороте. Я мог бы вручную создать вектор ограничений, но правда?


person a different ben    schedule 08.07.2014    source источник
comment
не используйте @ - попробуйте lim = as.vector (t (bbox (cp.proj))), но я не знаю, хорошо ли map.grid работает с пространственными объектами после построения, например, как map.grid знает что вам нужна сетка longlat, и в какой проекции находятся ваши данные?   -  person mdsumner    schedule 08.07.2014
comment
@mdsumner - спасибо. Понятия не имею, знает ли map.grid (), что делать! Мне сейчас сложно что-то делать. Любые предложения приветствуются.   -  person a different ben    schedule 08.07.2014
comment
@mdsumner -? map.grid () говорит: Draw a latitude/longitude grid on a projected map, но я полагаю, это говорит о карте, созданной функцией map.   -  person a different ben    schedule 08.07.2014
comment
Думаю, это правильно. Я уверен, что mapproj по умолчанию использует единичную сферу, так что вы собираетесь продвигать его в гору. См. «Линии сетки для работы с sp. (Это кошмар, который я знаю).   -  person mdsumner    schedule 08.07.2014
comment
Ага! rgdal имеет llgridlines, и я думаю, что это так! @mdsumner   -  person a different ben    schedule 08.07.2014
comment
Хорошо, кажется, что llgridlines, но маркировка исправлена. Здесь начинается дрожь и пот.   -  person a different ben    schedule 08.07.2014


Ответы (1)


РЕДАКТИРОВАТЬ: OP указывает rgdal::llgridlines, что является лучшим решением.

Вы используете контекст из sp/rgdal, который использует другую систему, чем mapproj/maps.

Попробуйте это (непроверено):

library(rgdal)
gl <- gridlines(myPoly)
cp.gl <- spTransform(gl, lcc)
plot(cp.gl, add = TRUE)

См. ?gridlines для получения дополнительной информации об использовании этого с ярлыками. Я считаю, что это хорошо работает, пока вы держитесь подальше от полярных карт.

person mdsumner    schedule 08.07.2014