| Finding pixel values of NETCDF image based on lat/lon [message #89066] |
Mon, 21 July 2014 20:46  |
Marta Yebra
Messages: 16
Registered: August 2012
|
Junior Member |
|
|
Hi all,
I am just introducing myself to NETCDF format and I have been able to read and access the data but Now I need to extract a time-serie of the values for a point location (lat/lon) and I do not know how to do it.
This is my piece of bad coding :):
filename='C:\UserData\yebram\GPP\2014516145028EnsembleGPP_GL .nc'
file_sites='C:\UserData\yebram\GPP\Sites_locations.csv'
;1. Open the file and assign it a file ID
NCDF_Id = ncdf_open(filename)
;When you are completely through with the file you should close it using the ncdf_close, fileID command.
;2. Find the number of file attributes and datasets (or variables). The information will be contained in the structure variable that we have named 'fileinq_struct',
;but you may give it any name you wish so long as you use the proper record names.
fileinq_struct=ncdf_inquire(NCDF_Id)
Ndims=fileinq_struct.Ndims ;The number of dimensions defined for this NetCDF file.
nvars = fileinq_struct.nvars ;The number of variables defined for this NetCDF file.
ngatts = fileinq_struct.ngatts ;The number of global attributes defined for this NetCDF file
RecDim=fileinq_struct.RecDim; The ID of the unlimited dimension, if there is one, for this NetCDF file. If there is no unlimited dimension, RecDim is set to -1.
; retrieve GPP data
; units = "kg m-2 s-1";
; _FillValue = -9999.0f; // float
; float gpp(time=360, lat=360, lon=720);
nameGPP='gpp'
NCDF_VARID_GPP = NCDF_VARID(NCDF_Id, nameGPP)
NCDF_VARGET, NCDF_Id, NCDF_VARID_GPP, GPP
GPP[where(GPP eq -9999)] = !VALUES.F_NAN
; retrieve long data
; units = "degrees_east";
nameLon='lon'
NCDF_VARID_lon = NCDF_VARID(NCDF_Id, nameLon)
NCDF_VARGET, NCDF_Id, NCDF_VARID_lon, Lon
Lon[where(Lon eq -9999)] = !VALUES.F_NAN
; retrieve lat data
; units = "degrees_north";
nameLat='lat'
NCDF_VARID_lat = NCDF_VARID(NCDF_Id, nameLat)
NCDF_VARGET, NCDF_Id, NCDF_VARID_lat, lat
lat[where(lat eq -9999)] = !VALUES.F_NAN
; retrieve time data
; units = "days since 1582-10-14 00:00:00";
; calendar = "standard";
nameTime='time'
NCDF_VARID_Time = NCDF_VARID(NCDF_Id, nameTime)
NCDF_VARGET, NCDF_Id, NCDF_VARID_Time, Time
Time[where(Time eq -9999)] = !VALUES.F_NAN
;open file with coordinates for teh sites
my_data=read_csvok(file_sites)
sites=my_data.Name_site
Lat_sites=my_data.Lat
Lon_sites=my_data.Lon
find_lon=where(Lon EQ Lon_sites, countlon)
find_lat=where(Lat EQ Lat_sites, countLat)
but this cannot find anything since my coordinate is for example 41.160000 and that value is not in the NETCDF.
How can I select the pixel that contains my point?
Your help is very much appreciate.
Kind regards,
Marta
|
|
|
|
|
|
|
|
| Re: Finding pixel values of NETCDF image based on lat/lon [message #89078 is a reply to message #89077] |
Tue, 22 July 2014 23:44  |
Fabzi
Messages: 305
Registered: July 2010
|
Senior Member |
|
|
On 23.07.2014 08:04, Marta Yebra wrote:
> Do I need to convert the NETCDF to GEOTIFF?
No!
In your example you write:
> sites = my_data.Name_site
> Lat_sites=my_data.Lat
> Lon_sites=my_data.Lon
> find_lon=where(Lon EQ Lon_sites, countlon)
> find_lat=where(Lat EQ Lat_sites, countLat)
There are two problems here:
- the EQ operator does not work when the arrays are of different length
- as you point out, where() cannot work since what you really need is
the coordinates of the nearest grid point
For this, several options are available, depending on what you want to
do (interpolate the gridded data to your point? Take the nearest grid
point) and depending on whether or not the lon/lat grid in your netcdf
file is regular or not...
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
