Построение интерполированных данных на карте

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

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

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

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

EDIT: Спасибо Пол за его очень полезный пост! Вот что у меня есть сейчас. Не удалось заставить 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")

Ответ 1

У меня есть ряд замечаний в вашем сообщении:

Использование кригинга

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

Проекция данных

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

Для выполнения расчетов в пакете 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 ( "+ proj и т.д." ) определяет целевую проекцию. Какой прогноз подходит, зависит от вашего географического положения и размера вашей области исследования.

Построение графика с помощью ggplot2

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

poly_df = fortify(poly_Spatial)

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

Создание сетки интерполяции

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

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

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