[R-sig-Geo] How do I read and plot .nc data in R?

[Michael Sumner]

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
>>
>
>

	[[alternative HTML version deleted]]