# Calculating polygon around point and extract gridded data

I'm trying to calculate `50km` circles around lat-long points and then extract gridded rain data (`1km x 1km`). However, for some reason my long-lat points are shown in the center of the U.S (picture related) and not where they are supposed to be. I have tried several projection transformations which more or less all deliver the same result. Please find attached a working example which shows the problem on a subsample.

``````library("daymetr")
library("ncdf4")
library("raster")
library("rasterVis")
library("sf")

coords <-structure(list(id = 2L, Longitude = -71.4836857,Latitude = 43.0202135), row.names = 1L, class = "data.frame")

# Downloads precipitation tile from Oak Ridge National Laboratory Severs into tempdir()
start = 2015,
end = 2015,
param = "prcp")

ncpath <- paste0(tempdir(),"/prcp_2015_11935.nc")
nc = nc_open(ncpath)

raster = raster(ncpath, band=4)

crs(raster)

dat_sf <- st_as_sf(coords, coords = c("Longitude", "Latitude"), crs = 4326) %>%
st_transform(3488) # transform to NAD83 California Albers, which uses metres as its unit.

# Buffer circle by 50km and transform to crs of raster object
dat_circles <- st_buffer(dat_sf, dist = 50000) %>%
st_transform("+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25 +a=60 +rf=6378137 +lat_2=45")

# Create Map of Raster Data and Circles
mapTheme <- rasterTheme(region=brewer.pal(8,"Greens"))
plt <- levelplot(raster, margin=F, par.settings=mapTheme)
plt + layer(sp.polygons(as_Spatial(dat_circles), fill = "red"))
``````

Your circles are okay, it seems the raster is in the wrong position.

Its bbox in its coordinates is:

``````> bbox(raster)
min     max
s1 2101250 2318250
s2  315500  566500
``````

and these are supposedly in this coordinate system:

``````> projection(raster)
 "+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25 +a=60 +rf=6378137 +lat_2=45"
``````

and if you find out where the extent of the raster is in lat-long you get:

``````> projectExtent(raster,"+init=epsg:4326")
class       : RasterLayer
dimensions  : 251, 217, 54467  (nrow, ncol, ncell)
resolution  : 0.05864168, 3.686555e-09  (x, y)
extent      : 69.54706, 82.27231, -90, -89.99999  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
``````

Just to confirm, lets see where the xmin, ymin point projects to:

``````> st_transform(
st_sfc(
st_point(c(2101250,315500)),
crs=projection(raster)),4326)

POINT (70.91257 -89.99999)
``````

Time to investigate the original data file.

Running `gdalinfo` shows it seems to have sensible corner coords near your circles:

``````Corner Coordinates:
Upper Left  ( 2101250.000,  566500.000) ( 71d46'39.16"W, 44d27'34.88"N)
Lower Left  ( 2101250.000,  315500.000) ( 72d50' 8.44"W, 42d13'10.67"N)
Upper Right ( 2318250.000,  566500.000) ( 69d 6'41.78"W, 43d44'56.38"N)
Lower Right ( 2318250.000,  315500.000) ( 70d15' 9.65"W, 41d32' 3.26"N)
Center      ( 2209750.000,  441000.000) ( 70d59'48.83"W, 42d59'58.50"N)
``````

but it also has a two matrices (not vectors) of lat and long coordinates...

``````Subdatasets:
SUBDATASET_1_NAME=NETCDF:"prcp_2015_11935.nc":lat
SUBDATASET_1_DESC=[251x217] latitude (64-bit floating-point)
SUBDATASET_2_NAME=NETCDF:"prcp_2015_11935.nc":lon
SUBDATASET_2_DESC=[251x217] longitude (64-bit floating-point)
``````

We can read these matrices in and generate the lat-long coordinates of each point of the grid.

``````> xy = data.frame(c(ncvar_get(nc,"lon")), c(ncvar_get(nc,"lat")))
> xy = st_as_sf(xy,coords=c(1:2))
> st_crs(xy)=4326
> plot(xy)
> axis(1)
``````

Looks okay. If we transform to what the raster says its projection is then:

``````> plot(st_transform(xy, projection(raster)))
``````

but if you look at that plot its not axis-aligned. This looks like a rotated grid.

Now it gets weird. Lets see where the first point in that lat-long grid goes to in raster projection:

``````> st_transform(xy[1,],projection(raster))
Simple feature collection with 1 feature and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 20.78145 ymin: 5.008069 xmax: 20.78145 ymax: 5.008069
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=25 +lat_2=45 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +a=60 +b=59.99999059286434 +units=m +no_defs
geometry
1 POINT (20.78145 5.008069)
``````

That's a really small number. This is supposed to be metres. Maybe we're near the origin? But no, even a point far away from that location in lat-long has a small number. That's odd. Something must be amiss with the raster projection...

``````> projection(raster)
 "+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25 +a=60 +rf=6378137 +lat_2=45"
``````

`+a=60`. What's that?

`````` +a Semimajor radius of the ellipsoid axis
``````

It thinks the semi-major radius of the earth is 60 METRES!??!

If you take that out and instead put in a spherical earth model with a sensible radius you get numbers that are at least reasonable, but probably not correct:

``````> p4 = "+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25  +r=6378137 +lat_2=45"
> projection(raster) = p4
> st_transform(xy[1,],projection(raster))
Simple feature collection with 1 feature and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 2212760 ymin: 532499.9 xmax: 2212760 ymax: 532499.9
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=25 +lat_2=45 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
geometry
1 POINT (2212760 532499.9)
``````

...and I wouldn't trust them, especially since the lat-long and projected grids looked rotated.

# Summary

• `raster` is not dealing with the CRS correctly
• CRS in NetCDF files are often weird
• NetCDF files will sometimes supply you with a lat and long matrix pair
• If it is a rotated coordinate system then you need to use the matrix pair
• The `stars` package will work with rotated coordinate system rasters

Also:

• I do not know where the `+a` parameter has come from and I can't find it in the `gdalinfo` metadata output. I suspect its confused it with one of the standard parallel values?
``````      Metadata:
lambert_conformal_conic#false_easting=0
lambert_conformal_conic#false_northing=0
lambert_conformal_conic#grid_mapping_name=lambert_conformal_conic
lambert_conformal_conic#inverse_flattening=298.257223563
lambert_conformal_conic#latitude_of_projection_origin=42.5
lambert_conformal_conic#longitude_of_central_meridian=-100
lambert_conformal_conic#semi_major_axis=6378137
lambert_conformal_conic#standard_parallel={25,60}
``````

It may be possible to construct a correct proj4 string from the metadata by hand.

Its also output correctly in the other `gdalinfo` metadata output:

``````PROJCS["unnamed",
GEOGCS["unknown",
DATUM["unknown",
SPHEROID["Spheroid",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",25],
PARAMETER["standard_parallel_2",60],
PARAMETER["latitude_of_origin",42.5],
PARAMETER["central_meridian",-100],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
``````

It may even be related to the warning I get when reading the raster:

``````> raster = raster(ncpath, band=4)
Warning message:
In cbind(m[i, ], vals) :
number of rows of result is not a multiple of vector length (arg 2)
``````

# Conclusion

Read the data in using `ncdf4` and create a curvilinear `stars` object from the lat-long matrix pair. Transform to another CRS from there.

# Update!

On my work machine, as opposed the my laptop where I did everything above, it seems to do the right thing:

``````> r = raster(ncpath, band=4)
> projection(r)
 "+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +a=6378137 +rf=298.257223563 +lat_1=25 +lat_2=60"
``````

that has the correct earth radius. When projected to lat-long, it seems to be in the right place approximately:

``````> plot(r)
> projectExtent(r,"+init=epsg:4326")
class       : RasterLayer
dimensions  : 251, 217, 54467  (nrow, ncol, ncell)
resolution  : 0.01716163, 0.01164861  (x, y)
extent      : -72.83568, -69.1116, 41.53589, 44.45969  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
``````

This is a slightly old R and these versions of the relevant packages:

``````rgdal_1.3-3
raster_2.6-7
R version 3.4.4 (2018-03-15)
``````

`sf` and `rgdal` report thus:

``````> library(sf)
Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3

> library(rgdal)
rgdal: version: 1.3-3, (SVN revision 759)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
Path to GDAL shared files: /usr/share/gdal/2.2
GDAL binary built with GEOS: TRUE
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: (autodetected)
• This PROJ4 from your update does the trick: `"+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +a=6378137 +rf=298.257223563 +lat_1=25 +lat_2=60"`. – cgx May 28 '19 at 14:04