Re: [R] Extracting point data using longitude and latitude from netcdf file using R

360 views
Skip to first unread message

Peter Tuju via R-help

unread,
Jan 9, 2016, 6:32:56 AM1/9/16
to r-h...@r-project.org
I have data file in netcdf with three dimensions (x, y, t) and I want to extract a variable RAINC and RAINNC
using longitude and latitude for a single point location with all the time, but no lucky. The syntax is as follows;;
setwd( "/run/media/tuju/0767090047/extract_wrf_txt_file" )
rm( list = ls() )                                             
library( ncdf4 )                                    
inp_file <- nc_open( "wrfout_d01_2016-01-07.nc" )time <- ncvar_get( inp_file, "Times" )                      # Reading the time variabledar_lon <- 39.2dar_lat <- -6.866667
RAINC <- ncvar_get( inp_file, varid = "RAINC", start  = c( dar_lon, dar_lat, 1 ), count = c( 1, 1, -1 ) )
RAINNC <- ncvar_get( inp_file, varid = "RAINNC", start  =  c( dar_lon, dar_lat, 1 ), count = c( 1, 1, -1 ) )
RAIN <- RAINC + RAINNCRAIN_TABLE <- cbind( time, RAIN )
write.table( RAIN_TABLE, "Dar_es_Salaam.txt", row.names = FALSE,
             col.names = c( "Valid Forecast Time",  "Rain (mm)", sep = "\t " )

# But no lucky with the red bolded syntax as I end up with the following error message> RAINC <- ncvar_get( inp_file, varid = "RAINC", start  = c( Lon[2], Lat[2], 1 ), count = c( 1, 1, -1 ) )
Error in Rsx_nc4_get_vara_double: NetCDF: Index exceeds dimension bound
Var: RAINC  Ndims: 3   Start: 0,4294967289,38 Count: 17,1,1
Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, addOffset,  :
  C function R_nc4_get_vara_double returned error
However when I cahnge the latitude to postive it works fine. Note latitudes in the file data ranges from -16.71505 to 7.787529 as shown below;
head(ncvar_get(inp_file, "XLAT"))
[1] -16.71505 -16.71505 -16.71505 -16.71505 -16.71505 -16.71505
> tail(ncvar_get(inp_file, "XLAT"))
[1] 7.787529 7.787529 7.787529 7.787529 7.787529 7.787529
## So, how can I get the syntax correct? Please help _____________
Peter  E. Tuju
Dar es Salaam
T A N Z A N I A
----------------------



[[alternative HTML version deleted]]

______________________________________________
R-h...@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Anthoni, Peter (IMK)

unread,
Jan 9, 2016, 7:16:36 AM1/9/16
to Peter Tuju, r-h...@r-project.org
Hi Peter,

the start in nc_varget requires a latitude and longitude index, not the latitude and longitude in double format.
So you need to figure out what index your latitude and longitude correspond to, which will depends on what data are in your netCDF.

it might have looked like that it worked for a positive latitude, but you got the data from the latitude index 6 or 7, depends on how the double was transformed into an integer.

best regards
Peter

Peter Tuju via R-help

unread,
Jan 9, 2016, 7:59:17 AM1/9/16
to Anthoni, Peter (IMK), r-h...@r-project.org
Thank you Mr. Anthon for your feedback. So, how can  I extract the data at the specified  latitude and longitude?Or, How can I get the latitude and longitude index which corresponds to the following coordinates from the netcdf file?
dar_lon <- 39.2
dar_lat <- -6.866667 _____________

Peter  E. Tuju
Dar es Salaam
T A N Z A N I A
----------------------

Ben Tupper

unread,
Jan 9, 2016, 8:02:32 AM1/9/16
to r-h...@r-project.org, Peter Tuju
Hi,

This post gives more details on how to transform your lat/lon values to row/column indices. The question and answer are specifically about the ncdf package, but the workflow is identical when using the ncfd4 package.

https://stat.ethz.ch/pipermail/r-help/2011-March/272641.html

Cheers,
Ben
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

Peter Tuju via R-help

unread,
Jan 9, 2016, 9:06:12 AM1/9/16
to Ben Tupper, r-h...@r-project.org
Thank you Be for the good guide, however no luck with the syntax used namely;
ix0 = wherenearest( lower_left_lon_lat[1],  lon )
ix1 = wherenearest( upper_right_lon_lat[1], lon )
iy0 = wherenearest( lower_left_lon_lat[2],  lat )
iy1 = wherenearest( upper_right_lon_lat[2], lat )
# I end up with this error,  "Error: could not find function "wherenearest"
Is there any other way I can get the index corresponding/or rearing to thelongitude and latitude of interests? _____________

Peter  E. Tuju
Dar es Salaam
T A N Z A N I A
----------------------

From: Ben Tupper <btu...@bigelow.org>
To: "r-h...@r-project.org" <r-h...@r-project.org>
Cc: Peter Tuju <pete...@ymail.com>
Sent: Saturday, January 9, 2016 4:01 PM
Subject: Re: [R] Extracting point data using longitude and latitude from netcdf file using R

Ben Tupper

unread,
Jan 9, 2016, 10:29:08 AM1/9/16
to r-h...@r-project.org
Hi,

Well, you have to extrapolate from that post that wherenearest is a function you must create on your own. Maybe something like this?

#' Find the index into a dataset that is 'closest' to a specified point.
#'
#' Adpated from https://stat.ethz.ch/pipermail/r-help/2011-March/272641.html
#' @param myPoint numeric, one point
#' @param allpoints numeric, one or more points
#' @return index into allPoints that myPoint falls closest to
wherenearest <- function(myPoint, allPoints){
d <- abs(allPoints-myPoint[1])
index <- which.min(d)
return( index )
}

I haven't tried the above.

You haven't provided much detail on what's in the file, but if it is a grid then perhaps you would find it easier to use the raster package. It has functions to read gridded data stored in ncdf files. It comes with a very good tutorial and makes working with gridded data a breeze. Using raster, you can bypass the nitty-gritty of getting data out of a ncdf file and just get to work on your data.

https://cran.r-project.org/web/packages/raster/index.html
https://cran.r-project.org/web/packages/raster/vignettes/Raster.pdf

Bon chance!
Ben

P.S. Do yourself (and everyone else on the list) a favor by making your email client use plain text rather than html or rich text when sending messages to the list. The html/rich text scrambles your code making it hard to read.
Reply all
Reply to author
Forward
0 new messages