[R-sig-Geo] How do I read and plot .nc data in R?
Michael Sumner
mdsumner at gmail.com
Fri Dec 19 10:04:56 CET 2014
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
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list