[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