Plotting Netcdf data. How to get right projection?

Question

Maybe it's very easy for someone to fix, but I'a a new one in R, so I really need help. I/m trying to plot data containing lat, lon and variable (nmf2). The grid isn't a regular one, so there are some difficulties with plotting. I've already read Plotting NetCDF file using lat and lon contained in variables and about the other method, but I still can't get the right image. Here is the code:

library(raster)
inputfile <- "F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC"

lat <- raster(inputfile, varname="PIERCEPOINT_NIGHT_LATITUDE")
lon <- raster(inputfile, varname="PIERCEPOINT_NIGHT_LONGITUDE")

plat <- rasterToPoints(lat)
plon <- rasterToPoints(lon)
lonlat <- cbind(plon[,3], plat[,3])

lonlat <- SpatialPoints(lonlat, proj4string = CRS("+proj=longlat +datum=WGS84"))
extent(lonlat)
#class       : Extent 
#xmin        : 0.008961686 
#xmax        : 359.983 
#ymin        : -84.95161 
#ymax        : 89.68419 

pr <- raster(inputfile, varname="NMF2_DISK")
proj4string(pr) <- "+proj=longlat +datum=WGS84"

pr
#class       : RasterLayer 
#dimensions  : 408, 13, 5304  (nrow, ncol, ncell)
#resolution  : 27.69031, 0.4280289  (x, y)
#extent      : 0.008961686, 359.983, -84.95161, 89.68419  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0  
#data source : C:\Users\Svetlana\Science\GUVI\R\SSUSI\F18-SSUSI_EDR-NIGHT- DISK_DD.20150107_SN.26920-00_DF.NC 
#names       : NMF2_DISK 
#zvar        : NMF2_DISK 

plot(pr)
r <- projectRaster(pr, crs=CRS("+proj=longlat +datum=WGS84"))
plot(r)

the results are like this^

The date is available here https://docs.google.com/file/d/0B6IqnlmRMSpcNFBXWWlha1JUUzQ/edit?usp=docslist_api

Answer 1

This is not really a raster in longitude latitude, it's just arrays of values (including longitude and latitude). You can deal with these explicitly like this:

f <- "F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC"

library(raster)
## treat these not as rasters, but as arrays of values 
## though raster() is extremely helpful in simplifying the use of ncdf/ncdf4
xyznames <- c("PIERCEPOINT_NIGHT_LONGITUDE", "PIERCEPOINT_NIGHT_LATITUDE", "NMF2_DISK")
lon <- raster(f, varname= xyznames[1])
lat <- raster(f, varname= xyznames[2])
dat <- raster(f, varname = xyznames[3])

## values drops the "raster" wrapper, just returns values in order as a vector
d <- cbind(values(lon), values(lat), values(dat))
## remap arbitrary values to [0,1] for a colour table
scl <- function(x) (x - min(x, na.rm = TRUE))/diff(range(x, na.rm = TRUE))
n <- 56
plot(d[,1:2], pch = 16, col = terrain.colors(n)[scl(d[,3]) * (n-1) + 1])

What happens next depends on what you need to do - it wouldn't make much sense to try to rasterize this in longlat since it's following a satellite swath line. Do you need to deal with the whole globe or are you interested only in a local area? Will you be merging multiple swaths? There are track-line following satellite projections but I am not sure that would help.

Answer 2

The projection you are setting (proj4string(pr) <- "+proj=longlat +datum=WGS84") is wrong. That should be the projection of the original data, not the target projection. You say that the data is NOT on a lat-lon grid, so you have to find out which CRS your data has. Looking at the .nc file, it's not contained there, so you have to look through the docs from the file source, I assume. If you can define the type of grid and its main parameters (e.g. projection origin), one can manually build the CRS easily. But without this info there is no way you can project, because you don't know the starting projection.

Right now your projectRaster() line is just projecting the data from lat-lon to lat-lon and doing strange things.

If you can find a Grads .ctl file describing how to plot this data, you can try to guess the projection parameters from that.

Answer 3

Unfortunately, the data has no projection, because it is satellite raw data captured during the flight around the earth. The first two subdatasets contain the lonlat values, but these are not understood by GDAL. So I tried to make a scatterplot out of the data:

I extracted the first three datasets to XYZ format:

gdal_translate -of XYZ NETCDF:"F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC":PIERCEPOINT_NIGHT_LATITUDE lat.xyz
gdal_translate -of XYZ NETCDF:"F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC":PIERCEPOINT_NIGHT_LONGITUDE lon.xyz
gdal_translate -of XYZ NETCDF:"F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC":NMF2_DISK nmf2.xyz

The files contain the local center coordinates of the raster field, WGS longitude, latitude and data values. I merged the columns manually in Libre Office:

"x" "y" "longitude" "latitude" "nmf2"
0.5 0.5 175.0537414551 44.6098899841 9.97E+036
1.5 0.5 173.9629058838 44.5663757324 9.97E+036
2.5 0.5 172.7252960205 44.5043334961 9.97E+036
3.5 0.5 171.3649902344 44.4205322266 9.97E+036
4.5 0.5 169.9060974121 44.3123893738 9.97E+036
5.5 0.5 168.3715362549 44.1780853271 9.97E+036
6.5 0.5 166.7823791504 44.0165939331 9.97E+036
7.5 0.5 165.1578521729 43.8276977539 9.97E+036
8.5 0.5 163.5154876709 43.6119537354 9.97E+036
9.5 0.5 161.8714599609 43.3707046509 9.97E+036
10.5 0.5 160.2409820557 43.106086731 9.97E+036
11.5 0.5 158.6387634277 42.8210525513 9.97E+036
12.5 0.5 157.0793457031 42.5194320679 9.97E+036
0.5 1.5 175.176864624 43.7407646179 926910.625

Note that 9.97E+036 is the NODATA value.

Finally I loaded the result as delimited text into QGIS. This is what the data looks like in EPSG:4326:

The lontiudes are between 0 and 360°, so I have reprojected the file to Pseudo Mercator EPSG:3857 and back to WGS84 to get the usual -180/+180° view.

To get the same with R, you might need some programming.