23

Вопрос по ggplot2, r – Построение интерполированных данных на карте

У меня есть данные обследования богатства видов, которые были получены на различных участках в Чесапикском заливе, США, и я хотел бы графически представить данные в виде «тепловой карты».

У меня есть массив данных координат широты / долготы и значений богатства, которые я преобразовал вSpatialPointsDataFrame и использовалautoKrige() функция из пакета automap для генерации интерполированных значений.

Во-первых, кто-нибудь может прокомментировать, правильно ли я реализуюautoKrige() функционировать?

Во-вторых, у меня проблемы с отображением данных и наложением карты региона. Альтернативно, я мог бы указать сетку интерполяции, чтобы отразить границы залива (как предложеноВот)? Любые мысли о том, как я могу это сделать и где я могу получить эту информацию? Поставка сетки вautoKrige() кажется достаточно легким


РЕДАКТИРОВАТЬ: Спасибо Полу за его супер полезный пост! Вот что у меня сейчас. Возникли проблемы при получении ggplot для приема как интерполированных данных, так и проекции карты:

require(rgdal)
require(automap)
#Generate lat/long coordinates and richness data
set.seed(6)
df=data.frame(
  lat=sample(seq(36.9,39.3,by=0.01),100,rep=T),
  long=sample(seq(-76.5,-76,by=0.01),100,rep=T),
  fd=runif(10,0,10))
initial.df=df

#Convert dataframe into SpatialPointsDataFrame
coordinates(df)=~long+lat

#Project latlong coordinates onto an ellipse
proj4string(df)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
#+proj = the type of projection (lat/long)
#+ellps and +datum = the irregularity in the ellipse represented by planet earth

#Transform the projection into Euclidean distances
project_df=spTransform(df, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")) #projInfo(type="proj")

#Perform the interpolation using kriging
kr=autoKrige(fd~1,project_df)
#Extract the output and convert to dataframe for easy plotting with ggplot2
kr.output=as.data.frame(kr$krige_output)
#Plot the output
#Load the map data for the Chesapeake Bay
cb=data.frame(map("state",xlim=range(initial.df$long),ylim=range(initial.df$lat),plot=F)[c("x","y")])

ggplot()+
  geom_tile(data=kr.output,aes(x=x1,y=x2,fill=var1.pred))+  
  geom_path(data=cb,aes(x=x,y=y))+
  coord_map(projection="mercator")
  • Error: User Rate Limit Exceededcoord_mapError: User Rate Limit ExceededrgdalError: User Rate Limit Exceeded

    от jslefche
  • Error: User Rate Limit ExceededCRS()Error: User Rate Limit ExceededspTransformError: User Rate Limit Exceeded

    от jslefche
  • Error: User Rate Limit Exceededutm proj4Error: User Rate Limit Exceededcoordinate systemError: User Rate Limit Exceeded

    от
  • Error: User Rate Limit Exceeded

    от
  • Error: User Rate Limit Exceeded

    от jslefche
  • пример кода дает мне пустой график с ggplot_0.9.3 в R 2.15 и ошибкой из графика в R 3.0.0 (нет подходящего метода для глубины)

    от David LeBauer
  • 47

    У меня есть ряд замечаний к вашему посту: Using kriging

    У меня есть ряд замечаний к вашему посту:

    Using kriging

    Я вижу, что вы используете геостатистику для построения своей тепловой карты. Можно также рассмотреть другие методы интерполяции, такие как сплайны (например, сплайны тонких пластин в пакете полей). Они делают меньше предположений о данных (например, стационарность), а также могут отлично визуализировать ваши данные. Сокращение количества допущений может помочь в случае, если вы отправите его в журнал, тогда у вас будет меньше объяснений для рецензентов. Вы также можете сравнить несколько методов интерполяции, если хотите, смотритеотчет, который я написал за несколько советов.

    Data projection

    Я вижу, что вы используете лат длинные координаты для кригинга. Эдзер Пебесма (авторgstat) отметил, что Не существует моделей вариограмм, подходящих для координат широт. Это связано с тем, что в последнее время расстояния не являются прямыми (т.е.евклиды), но над сферой, (т.е.Большие расстояния круга). Не существует ковариационных функций (или моделей вариограмм), которые действительны для сферических координат. Я рекомендую проецировать их, используяspTransform отrgdal пакет перед использованием автома

    Пакет rgdal используетбиблиотека проекций proj4 выполнить расчеты. Чтобы спроецировать ваши данные, вам сначала нужно определить их проекцию:

    proj4string(df) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
    

    Строка proj4 в правой части выражения выше определяет тип проекции (+proj) эллипсы, которые использовались (+ellps) и данные (+datum). Чтобы понять, что означают эти термины, вы должны представить Землю как картофель. Земля не идеально сферическая, это определяется эллипсами. Земля также не является идеальным эллипсоидом, но поверхность более неровная. Эта нерегулярность определяется датумом. Смотрите такжеэта статья в Википедии.

    Как только вы определили проекцию, вы можете использоватьspTransform:

    project_df = spTransform(df, CRS("+proj= etcetc"))
    

    где CRS (& quot; + proj и т. д.) определяет целевой прогноз. Какой прогноз подходит, зависит от вашего географического положения и размера вашей области обучения.

    Plotting with ggplot2

    Для добавления полигонов или полилиний в ggplot, пожалуйста, посмотрите документациюcoord_map, Это включает в себя пример использованияmaps пакет для построения границ страны. Если вам нужно загрузить, например, шейп-файлы для вашей учебной области, вы можете сделать это, используяrgdal, Помните этоggplot2 работает с data.frame, а неSpatialPolygons, Вы можете преобразоватьSpatialPolygons вdata.frame с помощью:

    poly_df = fortify(poly_Spatial)
    

    Смотрите такжеэта функция Я создал для построения пространственных сеток. Это работает непосредственно на SpatialGrids / Pixels. Обратите внимание, что вам нужно получить один или два дополнительных файла из этого хранилища (continuousToDiscrete).

    Creating interpolation grid

    Я создал automap для генерации выходной сетки, когда ничего не было указано. Это делается путем создания выпуклой оболочки вокруг точек данных и выборки 5000 точек внутри нее. Границы области прогнозирования и количество точек, выбранных в ней (и, следовательно, разрешение), совершенно произвольны. Для конкретного приложения форму области прогнозирования можно получить из многоугольника, используяspsample для выборки точек внутри многоугольника. Сколько точек для выборки и, следовательно, разрешение зависит от двух вещей:

    the kind of data you have, For example, if your data is very smooth, there is not much point in raising the resolution really high in comparison to this smoothness. Alternatively, if your data has many small scale strcutures, you need a high resolution. This is only possible ofcourse if you have the observations to support this high resolution. the density of data. If your data is more dense, you can raise the resolution.

    Если вы используете свою интерполированную карту для последующего анализа, важно правильно выбрать разрешение. Если вы используете карту исключительно для целей визуализации, это менее важно. Тем не менее, обратите внимание, что в обоих случаях слишком высокое разрешение может вводить в заблуждение относительно точности ваших прогнозов, и что слишком низкое разрешение не соответствует данным.