[GIS] Getting time-zones from latitude and longitude in raster using R

rtime

I have a raster with longitudes and latitudes for the US with meteorological data (for about 100k grid points). I would like to get the time zone for each point in order to transform the time to local time. Here is how my raster look like:

> raster
class       : SpatialPointsDataFrame 
features    : 95025 
extent      : -179.9979, 179.9918, 1, 85.33336  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
variables   : 1
names       : narr.a_221_20141002_1200_000.sub 
min values  :                                0 
max values  :                              100

This problem is somewhat similar to https://stackoverflow.com/questions/23414340/convert-to-local-time-zone-using-latitude-and-longitude, except that I have lots of points. The suggested solution there does not work as the Google API accepts only limited number of requests per day and I do not think that the brute force solution would be the most efficient.

I am an absolute beginner in geographical data, but familiar with R, so I need a solution which works in R.

Best Answer

As the number of points are way to much to request the timezone by lon/lat for each, the solution of @MickyT worked:

tz <- shapefile("us/tz_us.shp")
tz_polygon <- SpatialPolygons(tz@polygons, proj4string=tz@proj4string)
pts_over_tz <- over(raster, tz)

(Also note that my "raster" is actually a SpatialPointsDataFrame object, created by rasterToPoints(r, spatial = TRUE) where r is indeed a RasterLayer object.)

Related Solutions

[GIS] Extract values from raster using point data with buffer in R

If buffer is in meters the code works, but if I use degrees, it won't.

Indeed; the buffer cannot be specified in degrees. ?extract, buffer says: "If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters)."

Do you really want the buffer to be in degrees?

Is so, temporarily set the CRS to UTM, do the extract, and set it back.

cr <- crs(r)
crs(r) <- '+proj=utm +zone=10 +datum=WGS84'
e <- extract(r, xy, buffer=0.5)
crs(r) <- cr

Process Vector to Raster Faster with R – Speed Optimization Techniques

I tried to "parallelize" the function rasterize using the R package parallel in this way:

split the SpatialPolygonsDataFrame object in n parts
rasterize every part separately
merge all the parts into one raster

In my computer, the parallelized rasterize function took 2.75 times less than the no-parallelized rasterize function.

Note: the code below download a polygon shapefile (~26.2 MB) from the web. You can use any SpatialPolygonDataFrame object. This is only an example.

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Figure 2: Brazil Raster plot