Reprojecting GOES-16 data in R

Question

I'm trying to plot GOES-16 data with world boundaries. Data comes from GOES-16 files such as this.

Data comes with "GOES-R ABI fixed grid projection" reference system, which to my better knowledge is describe as this PROJ4 string: '+proj=geos +h=35786023.0 +a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs'

My first step is to create a raster layer:

library(raster)
rrqpe <- raster("OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc", varname='RRQPE')

I get some warnings about not being able to create a valid CRS, perhaps due software version.

Warning messages: 1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS: long_name=GOES-R ABI fixed grid projection perspective_point_height=35786023 semi_minor_axis=6356752.31414 sweep_angle_axis=x 2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS grid_mapping_name; false_easting; false_northing; scale_factor_at_projection_origin; scale_factor_at_central_meridian; standard_parallel; standard_parallel1; standard_parallel2; longitude_of_central_meridian; longitude_of_projection_origin; latitude_of_projection_origin; > straight_vertical_longitude_from_pole; longitude_of_prime_meridian; semi_major_axis; inverse_flattening; earth_radius; +proj; +x_0; +y_0; +k_0; +k_0; +lat_1; +lat_1; +lat_2; +lon_0; +lon_0; +lat_0; +lon_0; +lon_0; +a; +rf; +a

> rrqpe
class       : RasterLayer 
dimensions  : 5424, 5424, 29419776  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -0.5, 5423.5, -0.5, 5423.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc 
names       : ABI.L2..Rainfall.Rate...Quantitative.Precipitation.Estimate 
zvar        : RRQPE 

Then I set the CRS:

> rrqpe@crs <- sp::CRS('+proj=geos +h=35786023.0 +a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs')
> rrqpe
class       : RasterLayer 
dimensions  : 5424, 5424, 29419776  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -0.5, 5423.5, -0.5, 5423.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=geos +h=35786023.0 +a=6378137.0 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs 
data source : OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc 
names       : ABI.L2..Rainfall.Rate...Quantitative.Precipitation.Estimate 
zvar        : RRQPE 

Then I try to warp (reproject) the image to EPSG:4326

> rrqpe_wgs84 <- projectRaster(rrqpe, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84")

But something is wrong. The extent is tiny. I was expecting something like half world of data.

> rrqpe_wgs84
class       : RasterLayer 
dimensions  : 5436, 5436, 29550096  (nrow, ncol, ncell)
resolution  : 8.98e-06, 9.04e-06  (x, y)
extent      : -89.50005, -89.45123, -4.774882e-05, 0.04909369  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /tmp/RtmpfHzSHu/raster/r_tmp_2020-08-21_231812_20555_52886.grd 
names       : ABI.L2..Rainfall.Rate...Quantitative.Precipitation.Estimate 
values      : -49.68673, 49.87008  (min, max)

I want to plot the data (rrqpe_wgs84) on a world map like world["continent"]["geom"] from library spData. Can someone help me on this reprojection using R libraries? Although I'm not using here, I have rgdal installed. I'm trying to avoid using many file operations with gdalUtils on R 3.4.4.

Answer 1

This extent when you read the data in looks a bit ropey to me:

extent      : -0.5, 5423.5, -0.5, 5423.5  (xmin, xmax, ymin, ymax)

especially since the raster is 5424x5424 pixels. Looks like R isn't picking up the extent properly.

If I load it into QGIS and inspect it I see the extent:

Extent  -5434894.7009821739047766,-5434895.2182222492992878 :
         5434895.2182222492992878,5434894.7009821739047766

Which if I apply to your data (with your CRS as well) I do:

> e = extent(raster(xmn=-5434894.7009821739047766,ymn=-5434895.2182222492992878 , xmx=5434895.2182222492992878,ymx=5434894.7009821739047766))
> extent(rrqpe)=e

I can then transform:

> rrqpe_wgs84 <- projectRaster(rrqpe, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84")
Warning messages:
1: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
  309 projected point(s) not finite

with a bunch of warning messages because (I think) we are looking at an image of one side of the earth in a square raster, so there's no earth in the corners.

That has a reasonably extent in WGS84 lat-long:

> extent(rrqpe_wgs84)
class      : Extent 
xmin       : -165.8505 
xmax       : -13.49854 
ymin       : -67.26557 
ymax       : 67.36223 

The CRS as read by QGIS is different too:

> crs = "+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs +sweep=x"

If I assign that to your raster and transform:

> crs(rrqpe) = crs
> rrqpe_wgs84_crs <- projectRaster(rrqpe, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84")

It gives a different longitude range (because the +lon_0 parameter is different)

> extent(rrqpe_wgs84_crs)
class      : Extent 
xmin       : -151.3505 
xmax       : 1.001455 
ymin       : -67.26557 
ymax       : 67.36223 

The gdalinfo command (which I think you can run using the gdalUtils package from R) can tell you these things if you ask about the subdataset:

gdalinfo NETCDF:OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc:RRQPE
Driver: netCDF/Network Common Data Format
Files: OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc
Size is 5424, 5424
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["Spheroid",6378137,298.2572221]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",-75],
    PARAMETER["satellite_height",35786023],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    EXTENSION["PROJ4","+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"]]
Origin = (-5434894.700982173904777,-5434895.218222249299288)
Pixel Size = (2004.017315487541055,2004.017315487541055)

and lots more metadata including:

  RRQPE#resolution=y: 0.000056 rad x: 0.000056 rad
  RRQPE#scale_factor=0.00152602
  RRQPE#standard_name=rainfall_rate
  RRQPE#units=mm h-1
  RRQPE#valid_range={0,-6}

Anyway, the basic problem was R not getting the extent right.