Search code examples
rnetcdfr-rastermap-projections

How do I convert a projected netcdf file to raster for plotting in R

I have an netcdf file of lightning strikes that contains variables for both x, y and Lat Lon, The coordinate reference system for the xy data is albers_conical_equal_area and there are a number of additional variables false easting, false northing, etc. That go along with the projection. I'd like to convert the file to raster and I'm required to fill in a crs parameter. I'd like the raster displayed in the projected format but don't know how to properly format the crs entry.

Is this even possible (does raster only work with lat lon data?)

If so how do I convert the info in the netcdf file to work with the crs parameter in the raster command.

I've tried to get the attributes from the netcdf file but it's not clear how to use them.

I've used the code below to get the attributes from the crs variable but ...

grid_mapping_name <- ncatt_get(lightning, "crs", "grid_mapping_name")
standard_parallel <- ncatt_get(lightning, "crs", "standard_parallel")
longitude_of_central_meridian <- ncatt_get(lightning, "crs", 
"longitude_of_central_meridian")
latitude_of_projection_origin <- ncatt_get(lightning, "crs", 
"latitude_of_projection_origin")
false_easting <- ncatt_get(lightning, "crs", "false_easting")
false_northing <- ncatt_get(lightning, "crs", "false_northing")

I'm not sure how to apply this info to the raster command

strikes <- raster(t(lightning.array)
        , xmn=min(x)
        , xmx=max(x)
        , ymn=min(y)
        , ymx=max(y)
        , crs=CRS("missing code here ")
)

The file I'm using is a single day of lightning strikes across the US and I foresee importing a lot of these files in order to aggregate strikes in an area of interest for an entire year so I'd like to figure out this project issue on this file first.

Solution

  • Use the proj4string

    Try CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs").

    What you're looking for is the called the proj4string, which is the string of characters that encodes the coordinate reference system of your points.

    For example, CRS("+proj=longlat +datum=WGS84") is a CRS object in the sp package with a proj4string of "+proj=longlat +datum=WGS84". It's the familiar long/lat CRS that many data come in.

    Your data is in albers conical equal area, which I had to google the proj4sting for. Here's what I did for future reference. I googled "albers_conical_equal_area proj4string" which took me here. Then I grabbed the proj4string. That's it!