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

Question

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

У меня есть массив данных координат широты / долготы и значений богатства, которые я преобразовал в 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")

Answer 1

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

Используя кригинг

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

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

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

Печать с помощью ggplot2

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

poly_df = fortify(poly_Spatial)

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

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

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

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

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