4

I guess this will sound simple to you but I just cannot get the right plot within R. I have a netcdf4 file and would like to visualize it within R just like I would do in QGIS and ArcGIS.

sample data

I tried to use insights from similar existing posts but did not succeed. The image get squeezed as shown below:

library(raster)
mycrs <- CRS("+proj=lcc +lat_1=35 +lat_2=51 +lat_0=39 +lon_0=14 +k=0.684241 +units=m +datum=WGS84 +no_defs")
pr <- raster(file.nc4, varname="precipitationCal")
projection(pr) <- mycrs
pr

pr

class       : RasterLayer 
dimensions  : 3600, 1800, 6480000  (nrow, ncol, ncell)
resolution  : 0.1, 0.1  (x, y)
extent      : -90, 90, -180, 180  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=35 +lat_2=51 +lat_0=39 +lon_0=14 +k=0.684241 +units=m +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\3B-HHR.MS.MRG.3IMERG.20150405-S170000-E172959.1020.V03D.HDF5.dap.nc4 
names       : precipitationCal 
zvar        : precipitationCal 

plot(pr)

NB: my region of interest is lat=c(41,84) and lon=c(-142,-53)

enter image description here

file.nc4 attributes:

Grid.GridHeader: BinMethod=ARITHMETIC_MEAN;
Registration=CENTER;
LatitudeResolution=0.1;
LongitudeResolution=0.1;
NorthBoundingCoordinate=90;
SouthBoundingCoordinate=-90;
EastBoundingCoordinate=180;
WestBoundingCoordinate=-180;
Origin=SOUTHWEST;

EDIT: Using @mdsummer's suggestion below:

r <- '3B-HHR.MS.MRG.3IMERG.20150405-S170000-E172959.1020.V03D.HDF5.dap.nc4'
pr <- raster(r, varname="precipitationCal")
library(maptools);
data(wrld_simpl)
projection(wrld_simpl)
proj4string(pr)=proj4string(wrld_simpl)
plot(t(pr))
plot(wrld_simpl, add = TRUE)

I got this image but it is still not on the right latlon projection?.

enter image description here

3

The problem is you are assigning metadata to the object, and probably even overwriting the existing metadata. This happens here:

projection(pr) <- mycrs

You should first run this to see what raster thinks it is already:

projection(pr)

or just print the object out to get a fully summary:

pr

Going out on a limb, I think you should warp this raster to get to your target:

library(raster)
mycrs <- CRS("+proj=lcc +lat_1=35 +lat_2=51 +lat_0=39 +lon_0=14 +k=0.684241 +units=m +datum=WGS84 +no_defs")
pr <- raster(file.nc4, varname="precipitationCal")

ppr <- projectRaster(pr, crs = mycrs)  ## create a new object for good measure

BUT, warping is a destructive no-holds-barred re-modelling process, and so you may need to limit the extent of the target or crop the input or many other possible options. It also looks like your data is transpose y/x and so you probably need to fix that first - make sure you can get a sensible map, try this:

 pr <- raster(file.nc4, varname="precipitationCal")
 plot(t(pr))
 library(maptools);data(wrld_simpl)
 plot(wrld_simpl, add = TRUE)

(Untested since I can't access dropbox files right now).

If raster didn't already know what projection your data is in ("+proj=longlat +ellps=WGS84 etc. etc") then the projectRaster step would fail, and you would need to assign the right longlat PROJ.4 string before reprojecting.

  • @mdsummer thanks for providing feedback. You can access the files here disc.sci.gsfc.nasa.gov/SSW/#keywords=GPM_3IMERGHH 03; ignore the spatial bounding box tab. Date range can be 2015-06-01 to 2015-06-02. Just download one of the files. Choose precipitationCal only. – code123 Feb 11 '16 at 5:07

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.