Re: finding nearest point and 4 nearest points

Dorita,

For just the closest grid point, the function getind_latlon2d does
exactly that. See the example on the function page.

Do you really want the four closest points, or do you want the four
corners of the grid cell containing the given lat/lon? They are not
the same.

--Dave

On Thu, May 3, 2012 at 8:12 AM, Scott Capps <scapps@atmos.ucla.edu> wrote:
> I believe this would find the closest gridcells within 5km of stn_lat and
> stn_lon:
>
>   lon          = file2->XLONG(0,:,:)
>   lat          = file2->XLAT(0,:,:)
>   dist        = gc_latlon(stn_lat,stn_lon,lat,lon,2,4)
>   ndist       = dimsizes(dist)
>   dist1d      = ndtooned(dist)
>   tmp         = dimsizes(ind(dist1d.le.5.0))
>   if (tmp.lt.2) then
>      print("No gridcells found within 5km...")
>      exit
>   else
>      distidx  = ind_resolve(ind(dist1d.le.5.0),ndist) ; within 5km
>      ndistidx    = dimsizes(distidx)
>   end if
>   delete(tmp)
>
>
> Scott
>
> On 5/3/12 5:36 AM, Dorita Rostkier-Edelstein wrote:
>
> Hi folks,
>
> Does anyone have a script to find the closest and the 4 closest points to a
> given lat/lon, the nc file has the WRF XLAT, XLONG (pasted below).
>
> Thanks.
>
> Dorita
>
> netcdf wrf_static {
> dimensions:
>     Time = UNLIMITED ; // (1 currently)
>     south_north = 171 ;
>     west_east = 96 ;
> variables:
>     float HGT(Time, south_north, west_east) ;
>         HGT:FieldType = 104 ;
>         HGT:MemoryOrder = "XY " ;
>         HGT:description = "Terrain Height" ;
>         HGT:units = "m" ;
>         HGT:stagger = "" ;
>         HGT:coordinates = "XLONG XLAT" ;
>     float XLAND(Time, south_north, west_east) ;
>         XLAND:FieldType = 104 ;
>         XLAND:MemoryOrder = "XY " ;
>         XLAND:description = "LAND MASK (1 FOR LAND, 2 FOR WATER)" ;
>         XLAND:units = "" ;
>         XLAND:stagger = "" ;
>         XLAND:coordinates = "XLONG XLAT" ;
>     float XLAT(Time, south_north, west_east) ;
>         XLAT:FieldType = 104 ;
>         XLAT:MemoryOrder = "XY " ;
>         XLAT:description = "LATITUDE, SOUTH IS NEGATIVE" ;
>         XLAT:units = "degree_north" ;
>         XLAT:stagger = "" ;
>     float XLONG(Time, south_north, west_east) ;
>         XLONG:FieldType = 104 ;
>         XLONG:MemoryOrder = "XY " ;
>         XLONG:description = "LONGITUDE, WEST IS NEGATIVE" ;
>         XLONG:units = "degree_east" ;
>         XLONG:stagger = "" ;
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu May 3 15:34:29 2012