[Thread Prev][Thread Next][Index]
Re: how to output grid data into the arcview shape format file
Takaya,
Attached, please find the arc_ascii.jnl script which LAS uses to
generate an ASCII representation of data that can be ingested by ArcInfo
-- though not by ArcView. If anyone knows of a non-proprietary format
for creating gridded files for GIS systems, we'd love to know about it.
-- Jonathan Callahan
Namba Takaya wrote:
>
> Dir Sir
>
> Could somebody tell me how to output the grid data into the arc/view shape
> file by ferret?
>
> In the Live Access Server, there is a arc/view shape file format in the
> output format option.
> I wonder how the ferret(script) list the shape file out in the LAS server.
>
> Thank you in advance.
>
> Takaya Namba
>
> _________________________________________________________________
> $B%$%s%?!<%M%C%H$r$V$i$V$i%7%g%C%T%s%0$9$k$J$i(B MSN $B%7%g%C%T%s%0$X(B
> http://shopping.msn.co.jp/
!\cancel mode verify
! arc_ascii.jnl - create Arcinfo GIS output of regular grids
! Description: create Arcinfo GIS output of regular grids
! arguments: 1 2 3
! Usage: GO arc_ascii dset var file
! Hack to get quoted string with leading '/'
def symbol foo = "$1"
def symbol bar = "$3"
! check the arguments
query/ignore $foo%<Usage: GO std_gif_xy dset var file%
query/ignore $2%<Usage: GO std_gif_xy dset var file%
query/ignore $bar%<Usage: GO std_gif_xy dset var file%
can view
set view full
! initialize the data set
SET DATA ($foo)
let expr = missing($2,-1E+34)
let nodata = -1E+34
! THIS SCRIPT IS ONLY APPLICABLE TO THE XY PLANE
define symbol shape `expr,return=shape`
QUERY/IGNORE ($shape"|XY|<Region error: you can only obtain GIS output for XY maps")
set mode cal:days
set list/file="$3"
! subscript limits in XY plane
let istart `expr,return=istart`
let iend `expr,return=iend`
let jstart `expr,return=jstart`
let jend `expr,return=jend`
! size of grid
let Ni = iend-istart+1
let Nj = jend-jstart+1
! coordinate arrays
LET xcoords = x + 0*MISSING(expr[j=`jstart`],0)
LET ycoords = y + 0*MISSING(expr[i=`istart`],0)
! JC 990719 .. It looks like we don't need xmid,ymid ..
! coordinate bounds of grid
! Make sure West long is < 0, change both hi & lo if xlo is > 180
let slon = xcoords[i=`istart`]
IF `slon GT 180` THEN
let xlo = xcoords[i=`istart`] - 360
let xhi = xcoords[i=`iend`] - 360
ELSE
let xlo = xcoords[i=`istart`]
let xhi = xcoords[i=`iend`]
ENDIF
let ylo = ycoords[j=`jstart`]
let yhi = ycoords[j=`jend`]
! midpoint of grid
let xmid = (xlo+xhi) / 2
let ymid = (ylo+yhi) / 2
! size of grid cells
IF `Ni GT 2` THEN
let xcell_size = (xhi-xlo)/(Ni-1)
ELIF `Ni EQ 2` THEN
let xcell_size = (xhi-xlo)
ELSE
let xcell_size = 1
ENDIF
IF `Nj GT 2` THEN
let ycell_size = (yhi-ylo)/(Nj-1)
ELIF `Nj EQ 2` THEN
let ycell_size = (yhi-ylo)
ELSE
let ycell_size = 1
ENDIF
let cell_size = (xcell_size+ycell_size)/2
! SCRIPT IS ONLY APPLICABLE WHERE X and Y spacings are equal
IF `ABS(xcell_size-ycell_size)/(xcell_size+ycell_size) GT 1E-6` THEN
list/file/nohead/append/format="('WARNING: cell sizes in X and Y differ:',2(1PG15.7))" xcell_size, ycell_size
ENDIF
! make a heading as in this example
! ncols 180
! nrows 359
! xllcenter -89.54166
! yllcenter 0.4583334
! cellsize 1.
spawn touch ($bar)
spawn echo "ncols `Ni`" >> ($bar)
spawn echo "nrows `Nj`" >> ($bar)
list/file/nohead/append/format="('xllcorner',1PG15.7)" xlo
list/file/nohead/append/format="('yllcorner',1PG15.7)" ylo
list/file/nohead/append/format="('cellsize',1PG15.7)" cell_size
list/file/nohead/append/format="('nodata',1PG15.7)" nodata
can region/y
! reverse the j indexing
let jj = j[g=$2]
let j_rev = sortj(jj[j=`jstart`:`jend`]*-1)
let expr_rev = samplej(expr,j_rev)
! write out the data
list/file/nohead/append/format="(1PG15.7)" expr_rev
set mode/last verify
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement