How do I read and plot .nc data in R?
I had sent a second email that did not make it to the archives it seems,
here it is:
I think it works fine, I downloaded your file and:
library(raster)
f <- "yfin_20141219.nc"
r <- raster(f)
plot(r)
## check we get sensible orientation and scaling and offset and
interpretation etc. etc.
library(rworldmap)
data(countriesLow)
plot(countriesLow, add = TRUE)
r
: RasterLayer
dimensions : 541, 649, 351109 (nrow, ncol, ncell)
resolution : 0.08333333, 0.08333333 (x, y)
extent : 89.95833, 144.0417, -20.04167, 25.04167 (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=longlat +datum=WGS84
data source : C:\Users\mdsumner\Downloads\yfin_20141219.nc
names : yft_adu
z-value : 251164800
zvar : yft_adu
Note that in general, reading a raw array from file and applying this
assumptive georeferencing is fraught with danger:
?
y <- raster(data.adl, xmn=-20, xmx=25, ymn=90, ymx=144)
You should check assumptions about the orientation, order, units and
scaling of the coordinate arrays. The raster package is pretty smart about
telling you when danger lurks, but I would still be very careful. NetCDF
doesn't generally get used store the standard GIS affine transform, it
usually stores every dimension's coordinate and only removes that
redundancy in some cases for rectilinear grids. It's one of those murky
boundaries between different worlds.
HTH
Cheers, Mike.
On Fri Dec 19 2014 at 20:33:20 Eko Susilo <ekosusilo at live.com> wrote:
Hi Micheal,
I have tried do your scripts,
*y <- raster(fname, xmn=-20, xmx=25, ymn=90, ymx=144)*
* projection(y) <- CRS("+proj=longlat +datum=WGS84")*
the result is good, but not georeference image (i mean the data is not
right position), you can see my plot here
https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21163
Meanwhile when I read .nc directly and plot the data is look like
https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21164
<https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21163>
not georeferenced image but in the right position)
The .nc data available in this link
https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21165
On 19/12/2014 17:04, Michael Sumner wrote:
If you could easily convert it to raster, raster is the best place to
start.
Try
library(raster)
r = raster(fname)
plot(r)
Please report on the effect.
Netcdf is *so* general a format it really can come down to specifics. If
you can it is best to share a link to the or a representative file.
Cheers, Mike
On Fri, 19 Dec 2014 19:05 Eko Susilo <ekosusilo at live.com> wrote:
Could somebody help me to figure out my data
I have .nc data around Indonesia water (attached below)
I use R to open and plot it. This is my syntax
nc <- open.ncdf(choose.files('D:/R/yfin_pelikan/yfin_input', caption =
'Select netcdf Files .....'))
print (nc)
# Read the variable
adl <- nc$var[[1]]
varsize <- adl$varsize
ndims <- adl$ndims
nt <- varsize[ndims]
start <- rep(1,ndims) # begin with start=(1,1,1,...,1)
start[ndims] <- j # change to start=(1,1,1,...,i) to read timestep i
count <- varsize # begin w/count=(nx,ny,nz,...,nt), reads entire var
count[ndims] <- 1 # change to count=(nx,ny,nz,...,1) to read 1 tstep
data <- get.var.ncdf(nc, adl, start=start, count=count )
# plot the map
jet.colors <- colorRampPalette(c('#00007F', 'blue', '#007FFF',
'cyan','#7FFF7F', 'yellow', '#FF7F00', 'red', '#7F0000'))
plot.adl <- levelplot(data,
col.regions = jet.colors(255), at=seq(0, 1, 0.005),
yscale.components=yscale.raster.subticks,
xscale.components=xscale.raster.subticks,
margin=FALSE, ylab='latitude', xlab='longitude', main='test')
print(plot.adl + layer(sp.lines(indo.map, lwd=1, col='black')))
# The results is good bad not georreference so i can't overlay with the
base map (costaline)...
# So i try to convert it to raster
# Convert it to raster and re-plot
y <- raster(data.adl, xmn=-20, xmx=25, ymn=90, ymx=144)
projection(y) <- CRS("+proj=longlat +datum=WGS84")
plot.y <- levelplot(y,
col.regions = jet.colors(255), at=seq(0, 1, 0.005),
yscale.components=yscale.raster.subticks,
xscale.components=xscale.raster.subticks,
margin=FALSE, ylab='latitude', xlab='longitude', main='test')
print(plot.y + layer(sp.lines(indo.map, lwd=1, col='black')))
# The result is not good..
Attached file:
data (https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21165
<
https://www.researchgate.net/go.Deref.html?url=https%3A%2F%2Fonedrive.live.com%2Fredir%3Fresid%3D6FFDD661570C7D0A%2521165
)
result1 (https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21164 < https://www.researchgate.net/go.Deref.html?url=https%3A%2F%2Fonedrive.live.com%2Fredir%3Fresid%3D6FFDD661570C7D0A%2521164
)
result2 (https://onedrive.live.com/redir?resid=6FFDD661570C7D0A%21163 < https://www.researchgate.net/go.Deref.html?url=https%3A%2F%2Fonedrive.live.com%2Fredir%3Fresid%3D6FFDD661570C7D0A%2521163
)
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo