1

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.

GabrielMontenegro
  • 683
  • 1
  • 6
  • 21

1 Answers1

4

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

Val
  • 6,585
  • 5
  • 22
  • 52