Description
This example shows: how to use codes_grib_find_nearest and codes_get_element
Source code
grib_nearest.f90
! (C) Copyright 2005- ECMWF. ! ! This software is licensed under the terms of the Apache Licence Version 2.0 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. ! ! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by ! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. ! ! ! Description: how to use codes_grib_find_nearest and codes_get_element ! ! ! program find use eccodes implicit none integer :: npoints integer :: infile integer :: igrib, ios, i real(8), dimension(:), allocatable :: lats, lons real(8), dimension(:), allocatable :: nearest_lats, nearest_lons real(8), dimension(:), allocatable :: distances, values, lsm_values integer(kind=kindOfInt), dimension(:), allocatable :: indexes ! initialization open( unit=1, file="../../data/list_points",form="formatted",action="read") read(unit=1,fmt=*) npoints allocate(lats(npoints)) allocate(lons(npoints)) allocate(nearest_lats(npoints)) allocate(nearest_lons(npoints)) allocate(distances(npoints)) allocate(lsm_values(npoints)) allocate(values(npoints)) allocate(indexes(npoints)) do i=1,npoints read(unit=1,fmt=*, iostat=ios) lats(i), lons(i) if (ios /= 0) then npoints = i - 1 exit end if end do close(unit=1) call codes_open_file(infile, & '../../data/reduced_gaussian_lsm.grib1','r') ! A new grib message is loaded from file ! igrib is the grib id to be used in subsequent calls call codes_grib_new_from_file(infile,igrib) call codes_grib_find_nearest(igrib, .true., lats, lons, nearest_lats, nearest_lons,lsm_values, distances, indexes) call codes_release(igrib) call codes_close_file(infile) ! will apply it to another GRIB call codes_open_file(infile, & '../../data/reduced_gaussian_pressure_level.grib1','r') call codes_grib_new_from_file(infile,igrib) call codes_get_element(igrib,"values", indexes, values) call codes_release(igrib) call codes_close_file(infile) do i=1, npoints print*,lats(i), lons(i), nearest_lats(i), nearest_lons(i), distances(i), lsm_values(i), values(i) end do deallocate(lats) deallocate(lons) deallocate(nearest_lats) deallocate(nearest_lons) deallocate(distances) deallocate(lsm_values) deallocate(values) deallocate(indexes) end program find
grib_nearest.py
# # (C) Copyright 2005- ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. # import traceback import sys from eccodes import * INPUT = '../../data/reduced_gaussian_lsm.grib1' VERBOSE = 1 # verbose error reporting def example(): points = ((30, -20), (13, 234)) f = open(INPUT, 'rb') gid = codes_grib_new_from_file(f) for lat, lon in points: nearest = codes_grib_find_nearest(gid, lat, lon)[0] print(lat, lon) print(nearest.lat, nearest.lon, nearest.value, nearest.distance, nearest.index) four = codes_grib_find_nearest(gid, lat, lon, is_lsm=False, npoints=4) for i in range(len(four)): print("- %d -" % i) print(four[i]) print("-" * 100) codes_release(gid) f.close() def main(): try: example() except CodesInternalError as err: if VERBOSE: traceback.print_exc(file=sys.stderr) else: sys.stderr.write(err.msg + '\n') return 1 if __name__ == "__main__": sys.exit(main())
grib_nearest.c
/* * (C) Copyright 2005- ECMWF. * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. * * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ /* * C Implementation: nearest * * Description: how to get nearest point(s) */ #include <stdio.h> #include <stdlib.h> #include <string.h> #include "eccodes.h" static void usage(const char* prog) { printf("Usage: %s grib_file grib_file ...\n",prog); exit(1); } int main(int argc, char** argv) { int err = 0; long step=0; size_t nfiles; int i=0; codes_fieldset* set=NULL; codes_handle* h=NULL; char param[20]={0,}; size_t len=20; double lats[4]={0,}; double lons[4]={0,}; double values[4]={0,}; double distances[4]={0,}; int indexes[4]={0,}; char* order_by="param,step"; size_t size=4; double lat=-40,lon=15; int mode=0; int count; char** filenames; codes_nearest* nearest=NULL; if (argc < 2) usage(argv[0]); nfiles=argc-1; filenames=(char**)malloc(sizeof(char*)*nfiles); for (i=0;i<nfiles;i++) filenames[i]=(char*)strdup(argv[i+1]); set=codes_fieldset_new_from_files(0,filenames,nfiles,0,0,0,order_by,&err); CODES_CHECK(err,0); printf("\nordering by %s\n",order_by); printf("\n%d fields in the fieldset\n",codes_fieldset_count(set)); printf("n,step,param\n"); mode=CODES_NEAREST_SAME_GRID | CODES_NEAREST_SAME_POINT; count=1; while ((h=codes_fieldset_next_handle(set,&err))!=NULL) { CODES_CHECK(codes_get_long(h,"step",&step),0); len=20; CODES_CHECK(codes_get_string(h,"shortName",param,&len),0); printf("%d %ld %s ",count,step,param); if (!nearest) nearest=codes_grib_nearest_new(h,&err); CODES_CHECK(err,0); CODES_CHECK(codes_grib_nearest_find(nearest,h,lat,lon,mode,lats,lons,values,distances,indexes,&size),0); printf("\nIdx\tlat\tlon\tdist\tval\n"); for (i=0;i<4;i++) printf("%d\t%.2f\t%.2f\t%g\t%g\n", (int)indexes[i],lats[i],lons[i],distances[i],values[i]); printf("\n"); codes_handle_delete(h); count++; } if (nearest) codes_grib_nearest_delete(nearest); if (set) codes_fieldset_delete(set); return 0; }