0

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.

0

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.

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.