Создайте квадратные сетки и экспортируйте их как шейп-файл или таблицу - PullRequest
0 голосов
/ 03 апреля 2019

Я хотел бы использовать Google Earth Engine для извлечения данных для определенных стран. Мне нужны данные в виде квадратных сеток, поэтому я хотел бы создать эти квадратные сетки для определенной страны, добавить их в шейп-файл и затем импортировать шейп-файл в Earth Engine. Я уже нашел некоторый код для создания квадратных сеток ( Создание сетки внутри шейп-файла ), но теперь у меня есть две проблемы.

Во-первых, мне нужно экспортировать квадратные сетки, чтобы я мог импортировать их в Earth Engine. Я очень открыт для альтернатив шейп-файлу.

Во-вторых, последующий код работает для некоторых стран (например, Франции), но не для других (например, Таиланда).

library(raster)
shp = getData(country = "FRA", level = 0)

shp = spTransform(shp, CRSobj = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(shp)

cs = c(10000, 10000)
grdpts = makegrid(shp, cellsize = cs)

spgrd = SpatialPoints(grdpts, proj4string = CRS(proj4string(shp)))

spgrdWithin = SpatialPixels(spgrd[shp,])
plot(spgrdWithin, add = T)

Замена «FRA» на «THA» в строке 2 приводит к ошибке в spTransform .

1 Ответ

0 голосов
/ 04 апреля 2019

Сбой, потому что вы используете зону 32 utm. Вам нужно использовать зону исходя из долготы страны.Вы можете видеть их здесь

Вы можете автоматизировать поиск зоны с помощью ceiling((longitude+180)/6)

library(raster)
s <- getData(country = "FRA", level = 0)

Получить центроид.В этом случае вы можете сделать

centr <- coordinates(s)

Если есть несколько полигонов, вы можете сделать что-то вроде этого

centr <- apply(coordinates(s), 2, mean)

Вычислить UTM-зону.(обратите внимание, что у вас было 32 для Франции, что нехорошо)

zone <- ceiling((centr[1] + 180)/6)
zone
#[1] 31

А затем используйте это вот так

crs <- paste0("+proj=utm +datum=WGS84 +unit=m +zone=", zone)
st <- spTransform(s, crs)

Для Таиланда вы получите

s <- getData(country = "THA", level = 0)
centr <- apply(coordinates(s), 2, mean)
zone <- ceiling((centr[1] + 180)/6)
zone
#[1] 47

Однако это , а не подход, который будет работать для всех стран.Зоны UTM имеют ширину 6 градусов, и многие страны охватывают несколько зон (Россия берет пирог с 28 зонами).Поэтому, в зависимости от ваших целей, вы можете использовать другую систему координат (crs).

После этого альтернативный способ получения квадратных многоугольников состоит в создании слоя RasterLay со степенью s и разрешением по выбору.Но я сомневаюсь, что это лучший способ получить данные из GEE.Я бы предложил вместо этого загрузить схему страны.

r <- raster(st, res=10000)
r <- rasterize(st, r, 1)
x <- as(r, "SpatialPolygons")

# write to file
shapefile(x, "test.shp")

# view
plot(x)

enter image description here

...