I have several locations in longitude and latitude coordinates and I need to get an estimate of solar radiation for each location.
I found data from NASA that has exactly this information in a 1 degree x 1 degree resolution.
nasafile <- 'http://eosweb.larc.nasa.gov/sse/global/text/global_radiation'
nasa <- read.table(file=nasafile, skip=13, header=TRUE)
head(nasa)
Lat Lon Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Ann
1 -90 -180 9.63 5.28 0.75 0 0 0 0 0 0.1 3.24 8.28 10.97 3.19
2 -90 -179 9.63 5.28 0.75 0 0 0 0 0 0.1 3.24 8.28 10.97 3.19
3 -90 -178 9.63 5.28 0.75 0 0 0 0 0 0.1 3.24 8.28 10.97 3.19
4 -90 -177 9.63 5.28 0.75 0 0 0 0 0 0.1 3.24 8.28 10.97 3.19
5 -90 -176 9.63 5.28 0.75 0 0 0 0 0 0.1 3.24 8.28 10.97 3.19
6 -90 -175 9.63 5.28 0.75 0 0 0 0 0 0.1 3.24 8.28 10.97 3.19
As you can see, it has monthly data as well as an mean annual average for a paticular lat/lon coordinate in by 1 degrees.
If I convert this to a SpatialPixelsDataFrame
object:
proj <- CRS('+proj=latlon +ellps=WGS84')
coords <- nasa[,2:1]
data <- nasa[,3:15]
nasaSP <- SpatialPixelsDataFrame(points=coords, data=data, proj4string=proj)
I would like to know how to extract the annual value for a particular latitude and longitude coordinate.
I have no experience in geographic data analysis and I have been following this post: https://www.r-bloggers.com/maps-of-solar-radiation/ , but the point of it was to get a map and I just need to extract particular values and don't know how to do it.
I would convert the SpatialPixelsDataFrame
to a raster object in order to use the functionality of the raster
package. Particularly for extracting values, it's easier if it's a raster
object since you can extract values at any given location.
In this case, I'll convert it to a brick, which is a raster with multiple layers (since we have multiple variables).
library(raster)
# convert to brick
rbck <- brick(nasaSP)
# output
> rbck
class : RasterBrick
dimensions : 180, 360, 64800, 13 (nrow, ncol, ncell, nlayers)
resolution : 1, 1 (x, y)
extent : -180.5, 179.5, -90.5, 89.5 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, Ann
min values : 0.00, 0.00, 0.17, 0.00, 0.00, 0.00, 0.00, 0.00, 0.10, 0.00, 0.00, 0.00, 1.51
max values : 9.97, 8.20, 7.54, 7.79, 8.31, 9.05, 8.49, 7.89, 7.38, 8.28, 9.31, 11.19, 6.98
If you only want to look at one variable (you mentioned the annual value), you can easily subset it:
solAnn <- rbck[['Ann']]
# plot
plot(solAnn)
Now if you want to extract values, you can do it in multiple ways.
Suppose we have a matrix of coordinates:
extcoords <- matrix(c(1:10,1:10),ncol=2)
You can extract values either with the cell indices of the raster, which you can get with the function cellfromXY
:
vals <- solAnn[cellFromXY(solAnn,extcoords)]
> vals
[1] 5.49 5.43 5.28 5.06 4.53 4.81 5.12 5.37 5.52 5.79
Or you can convert the coordinates to spatial points and then use extract
to get your values:
pj <- crs('+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0')
pts <- SpatialPoints(coords = extcoords,proj4string = pj)
vals <- extract(solAnn,pts)
# Output
> vals
[1] 5.49 5.43 5.28 5.06 4.53 4.81 5.12 5.37 5.52 5.79
Using this method you can also add the points to your map with points()
.
If you're interested in all values, you can also extract directly from the raster brick before you subset:
vals <- extract(rbck,pts)
# Output
> vals
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Ann
[1,] 5.48 5.97 5.87 5.77 5.13 5.15 5.38 5.48 5.70 5.38 5.40 5.19 5.49
[2,] 5.63 6.12 5.95 5.81 4.96 4.76 4.98 5.27 5.40 5.29 5.54 5.39 5.43
[3,] 5.86 6.09 5.83 5.54 4.72 4.15 4.49 5.07 5.07 5.27 5.63 5.60 5.28
[4,] 5.84 6.01 5.67 5.40 4.62 3.78 3.93 4.66 4.53 5.00 5.58 5.65 5.06
[5,] 5.33 5.33 5.19 5.05 4.55 3.75 3.30 3.80 3.65 4.26 4.97 5.13 4.53
[6,] 5.56 5.60 5.39 5.16 4.86 4.44 3.94 3.82 4.04 4.50 5.06 5.32 4.81
[7,] 5.82 5.89 5.78 5.44 5.11 4.68 4.40 4.19 4.41 4.77 5.39 5.62 5.12
[8,] 5.90 6.15 6.15 5.75 5.35 4.88 4.48 4.28 4.58 5.13 5.87 5.87 5.37
[9,] 5.78 6.26 6.30 5.83 5.49 5.15 4.70 4.40 4.96 5.52 6.04 5.78 5.52
[10,] 5.69 6.15 6.40 6.24 6.08 5.70 5.35 5.08 5.46 5.84 5.85 5.59 5.79
HTH