2

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.

example

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()
download_daymet_tiles(location = c(43.59510,  -71.93360, 42.44533,  -71.0337),
                      start = 2015,
                      end = 2015,
                      param = "prcp")

## Read netCDF precipation data
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"))
1

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)
[1] "+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)
[1] "+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)
[1] "+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)
 Linking to sp version: 1.3-1 

I can compare with my laptop later today.

  • 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 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.