0

I would like to use the Google Earth Engine in order to extract data for certain countries. I need the data in the form of square grids, so I would like to create those square grids for a certain country, add them to the shapefile and then import the shapefile into the Earth Engine. I already found some code to create the square grids (Create a grid inside a shapefile), but now I have two problems.

First, I need to export the square grids so that I can import them to the Earth Engine. I'm very open to alternatives to the shapefile.

Second, the subsequent code works for some countries (like France), but not for others (like Thailand).

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)

Replacing "FRA" by "THA" in line 2 leads to an error in spTransform.

tho_mi
  • 698
  • 3
  • 9
  • 16

1 Answers1

0

That fails because you are using utm zone 32. You need to use the zone based on the longitude of the country. You can see them here

You can automate finding a zone with ceiling((longitude+180)/6) 

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

Get the centroid. In this case you can do

centr <- coordinates(s)

If there are multiple polygons, you can do something like this

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

Compute the UTM zone. (note that you had 32 for France, which is not good) 

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

And then use it like this

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

For Thailand you would get

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

However, this is not an approach that would work for all countries. UTM zones are 6 degrees wide, and many countries span multiple zones (Russia takes the cake with 28 zones). So depending on your goals, you may want to use another coordinate reference system (crs).

After that, an alternative way to get square polygons is to create a RasterLayer with the extent of s, and a resolution of choice. But I doubt that this is the best way to get data out of GEE. I would suggest uploading the country outline instead. 

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

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thanks for your answer! What approach would you suggest? Upload the outlines and then create the square grid in GEE? Or upload the outlines, download the image and do the rest in R? – tho_mi Apr 04 '19 at 20:03
  • I do not know what you want to accomplish in GEE, but either way, that is a separate question – Robert Hijmans Apr 05 '19 at 06:17
  • Then I will open a new one. Thanks for your help! – tho_mi Apr 05 '19 at 16:03
  • Additional question @Robert Hijmans. If I have a shapefile with a different projection ("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") the line "x <- as(r, "SpatialPolygons")" fails. Why? Translated error would be something like "no method or default for coercing "NULL" to "SpatialPolygons"). – tho_mi Apr 08 '19 at 10:02
  • I rephrase my question. I assume that it fails because resolution isn't in meters if the projection is longlat. What would be the way to go if I want to use meters, but not utm projection (in order to avoid the problem with zones)? – tho_mi Apr 08 '19 at 10:16