[ncl-talk] coordinates in raster file

Rick Brownrigg brownrig at ucar.edu
Mon Jun 3 19:22:48 MDT 2019


Hi,

As far as I can tell from visual inspection of the program logic and
printing values from the constructed arrays, this looks to be correct. Are
you seeing something that makes you doubt it?  I'm suspect about the need
to make 2D coordinate arrays, as, i) the grid is clearly rectilinear, and
ii) the docs for gsn_coordinates state it will take either 2D or 1D values:

    http://ncl.ucar.edu/Document/Graphics/Interfaces/gsn_coordinates.shtml

Hope that helps...
Rick


On Mon, Jun 3, 2019 at 5:24 PM Vanúcia Schumacher <
vanucia-schumacher at hotmail.com> wrote:

> Hi users,
>
> Could anyone check if I am reading and associating correctly the lat and
> lon coordinates of this raster file (attachment) with NCL?
>
> Thanks
>
> Part of the script:
>
>       files  = systemfunc("ls test.asc")
>        hdr = readAsciiHead(files, 6)
>      ncols = stringtoint  ( str_get_field(hdr(0), 2, " ") )
>      nrows = stringtoint  ( str_get_field(hdr(1), 2, " ") )
>      lonLL = stringtofloat( str_get_field(hdr(2), 2, " ") )
>      latLL = stringtofloat( str_get_field(hdr(3), 2, " ") )
>   deltaLon = stringtofloat( str_get_field(hdr(4), 2, " ") )
>   deltaLat = stringtofloat( str_get_field(hdr(4), 2, " ") )
> missingVal = stringtofloat( str_get_field(hdr(5), 2, " ") )
>
> data     = readAsciiTable(files,ncols,"float",6)
>                    data at _FillValue = missingVal
>
> ;----- construct 1D coordinates....
> lons             = ispan(0,ncols-1,1) * deltaLon + lonLL
> lats             = ispan(0,nrows-1,1) * deltaLat + latLL
> lats             = lats(::-1)
>
> ;---- gsn_coordinates wants 2D coordinate arrays
> data at lon2d       = conform_dims((/ nrows, ncols/), lons, 1)
> data at lat2d       = conform_dims((/ nrows, ncols/), lats, 0)
>
>   data!0                 = "lat"
>   data&lat              = lats
>   data&lat at units  = "degrees-north"
>   data!1                 = "lon"
>   data&lon            = lons
>   data&lon at units = "degrees-east"
>
> ;---- Open shapefile and read lat/lon values.
>   ...
>    data_mask = shapefile_mask_data(data,shp_filename,True)
>                         copy_VarMeta(data, data_mask)
>   data_point := (/data_mask({-34.64},{-70.32})/)
> ...
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190603/041ff847/attachment.html>


More information about the ncl-talk mailing list