CODEWITHSUNDEEP
rgisr-raster
Seb
Seb

Reputation: 5497

Extracting data from a NetCDF file with package raster in R and finding the proper projection

I downloaded a NetCDF-File here and read it via:

library(raster)
b <- brick("precipitationintensity14rcp45modelsensemblemedianannual.nc")

Then I checked the proj4string in order to see the projection and how to provide coordinates:

> proj4string(b)
[1] "+proj=lcc +lat_0=47.5 +lon_0=13.3299999237061 +lat_1=49 +lat_2=46 +x_0=400000 +y_0=400000 +datum=WGS84 +units=m +no_defs"

In order to extract data at a certain point I use long/lat coordinates in WGS84 format:

coords<-data.frame(lon=47.4, lat=11.5)
coordinates(coords)<-c("lon","lat")

Then I try to extract data by:

> extract(b$X1981.07.15,coords)
     [,1]
[1,]   NA

This leaves me somehow baffeled. Since I expected a value at this given point.

In order to debug my code I try to compare coordinates of the raster object with mine:

> coordinates(b$X1981.07.15)
               x      y
     [1,] 111500 571500
     [2,] 112500 571500
     [3,] 113500 571500

Now I'm confused since the coordinates provided do not match the proj4string format (WGES84). Could you please help me with:

  1. Providing propper coordinates for extracting data
  2. Explaining me why the coordinates in the raster object are not as I expected

Upvotes: 0

Views: 637

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47026

You need to transform the points to the coordinate reference system of the raster (not the other way around!). Below I show how to do that with "terra", the replacement of the "raster" package.

library(terra)
b <- rast("precipitationintensity14rcp45modelsensemblemedianannual.nc")

coords <- data.frame(lon=47.4, lat=11.5)
pts <- vect(coords, crs="+proj=longlat")

#crs(b) = "+proj=lcc +lat_0=47.5 +lon_0=13.3299999237061 ...
lccpts <- project(pts, crs(b))

extract(b$X1981.07.15, lccpts)

Upvotes: 1

Related Questions