Выравнивание большого растра и пространственных точек в R - PullRequest
Новая
       42

Выравнивание большого растра и пространственных точек в R

0 голосов
/ 02 ноября 2019

У меня есть два набора данных:

  1. A растровый файл , содержащий попиксельные классификации растительности по всей Австралии.

Описывается, что файл имеет проекцию EPSG: 3577 https://epsg.io/3577. Все дополнительные метаданные вместе с самими данными можно найти в http://www.agriculture.gov.au/abares/forestsaustralia/forest-data-maps-and-tools/spatial-data/forest-cover

A csv файл , содержащий данные столбца широты, столбца долготы и других переменных.

Поскольку это также большой набор данных, я предоставлю только заголовокиз этого: 

firms <- structure(list(latitude = c(-20.5897, -22.4119, -21.132, -23.9083, 
-23.5908, -23.689), longitude = c(147.6421, 148.8484, 148.1897, 
147.2966, 150.1676, 150.0994), brightness = c(320.7, 316.7, 320.5, 
320.1, 315.8, 314.3), scan = c(1.6, 2.1, 1.8, 1.7, 2.7, 2.7), 
    track = c(1.3, 1.4, 1.3, 1.3, 1.6, 1.6), acq_date = c("2017-01-01", 
    "2017-01-01", "2017-01-01", "2017-01-01", "2017-01-01", "2017-01-01"
    ), acq_time = c(47L, 47L, 47L, 47L, 47L, 47L), satellite = c("Terra", 
    "Terra", "Terra", "Terra", "Terra", "Terra"), instrument = c("MODIS", 
    "MODIS", "MODIS", "MODIS", "MODIS", "MODIS"), confidence = c(34L, 
    26L, 36L, 53L, 33L, 22L), version = c(6.2, 6.2, 6.2, 6.2, 
    6.2, 6.2), bright_t31 = c(299.6, 295.3, 299.7, 296.6, 291.7, 
    289.3), frp = c(19, 25.8, 22.9, 17.6, 35.4, 30), daynight = c("D", 
    "D", "D", "D", "D", "D"), type = c(0L, 0L, 0L, 0L, 0L, 0L
    )), row.names = c(NA, -6L), class = c("data.table", "data.frame"
), .internal.selfref = <pointer: 0x10200e2e0>)

Моя цель проста:

Построить растровые данные, а затем построить каждую строку из файла CSV в качестве точки на вершине растра.

Я доволен plot(), levelplot() или ggplot(), хотя последний кажется глупым, поскольку требуется слишком много времени для преобразования этого массивного растра в фрейм данных.

Пока мой код выглядит так: 

library(rgdal)
library(raster)
library(rasterVis)
library(sp)
  1. Загрузить в растровый файл 
test <- raster('w001000.adf')
Затем создайте координаты для data.frame и сделайте его форматом SpatialPoints 
coordinates(firms) <- ~longitude + latitude
firms_pts <- SpatialPoints(coords=coordinates(firms), proj4string = CRS("+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
Преобразовать данные точек в тот же crs, что и растр 
firms_lyr <- spTransform(firms_pts, crs(test))

Когда я пытаюсь построить данные в этой точке, я не вижу ни одной из точек.

Итак, я использую одинаковый экстент в обоих наборах данных, а затем строю график примерно друг на друге. 
extent(test) <- extent(firms)
сюжет 
plot(test)
points(firms_lyr, pch=16, cex=0.1)

СЕЙЧАС ... досадно, что наборы данных не выровнены должным образом пространственно this

Очевидно, что это что-тоделать с проекцией.

Я почти уверен, что делаю что-то не так в шагах 3 и / или 4

1 Ответ

1 голос
/ 03 ноября 2019

Это работает? Я не уверен, что это точки на Дальнем Востоке или нет.

data 

firms <- structure(list(latitude = c(-20.5897, -22.4119, -21.132, -23.9083,  -23.5908, -23.689), 
                        longitude = c(147.6421, 148.8484, 148.1897, 147.2966, 150.1676, 150.0994)), 
                   .Names = c("longitude", "latitude"), 
                   class = "data.frame", row.names = c(NA, -6L))

xy <- firms[,c(2,1)]
firms_pts <- SpatialPointsDataFrame(coords = xy, data = firms,
                       proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

Перепроецирование шейп-файла с использованием ADF 

library(raster)
setwd(paste(dirname(rstudioapi::getActiveDocumentContext()$path), "/aus_for18", sep=""))
test <- raster("aus_for18.adf")
firms_lyr <- spTransform(firms_pts, crs(test))

Plotting 

plot(test)
points(firms_lyr, pch=16, cex=1) #Queensland points

enter image description here

...