R
and Fortran
will read an ArcInfo binary grid directly, with essentially no extra processing, so skip the middleman: you don't need QGIS or anything else.
These datasets consist of two files. One is an ASCII header file with .hdr
extension, formatted to look like this:
ncols 1133
nrows 1415
xllcorner 686280
yllcorner 4179990
cellsize 30
NODATA_value -9999
byteorder LSBFIRST
The meanings should be evident.
The other file is an array of four-byte values (either floats or integers, depending on the type of the grid) arranged row by row, from top down.
Here is a quick and dirty example of how it can be read and processed in R
. (It's quick and dirty because it does not check for deviations from the format or inconsistencies, has no error checking, and has not been tested on all possible systems.) To demonstrate that the data manipulation is correct, I plot the data and show their coordinates.
Automation is trivial, of course: encapsulate this within a function and enclose it within an appropriate iterator (*apply
, for
, while
, or whatever in R
; FOR
, DO
, WHILE
, or whatever in Fortran).
R
's raster
package supports subsetting of grids by general polygons. For rectangular subsetting, that's a basic part of R
's array indexing capabilities.
#
# Read the header file.
#
sFile <- "G:/USGS/DEM/7_5min/VA/albem_s" # Root dataset name, without extensions
h <- read.table(paste(sFile, ".hdr", sep=""),
col.names=c("field", "value"), as.is=TRUE)
#
# Extract the header values.
#
header <- type.convert(h[1:6,"value"]); names(header) <- h[1:6,"field"]
n.rows <- header["nrows"]
n.cols <- header["ncols"]
nodata <- header["NODATA_value"]
x0 <- header["xllcorner"]
y0 <- header["yllcorner"]
cellsize <- header["cellsize"]
byteorder <- ifelse(h[7, "value"]=="LSBFIRST", "little", "big")
#
# Read the data.
#
x <- readBin(paste(sFile, ".flt", sep=""),
"numeric", n=n.rows*n.cols, size=4, endian=byteorder)
x[x==nodata] <- NA # Mark the NoData values
x <- t(matrix(x, nrow=1133)) # (R stores data by columns, not rows)
#
# Display the grid.
#
library(raster)
plot(raster(x, xmn=x0, xmx=x0+n.cols*cellsize, ymn=y0, ymx=y0+n.rows*cellsize))
This image looks correct. The total time to read and process the dataset (before displaying it) was 1/8 second.
One way to extract the points from adf is to open this adf file in QGIS and export it to a .xyz format and then use this file to match ur point to the point in the .xyz file
other way is to convert this adf file to a tif file using QGIS and then import this file to POSTGIS database using raster2pgsql command line tool.
this is how u import raster data (tif) to postgis
Install postgresql with postgis plugin
open command line and navigate to bin folder of postgresql normaly is it located here C:\Program Files\PostgreSQL\9.3\bin
STEP 1:
raster2pgsql -F -I -M -C "PATH of tif file with extension" public.your_tablename > your_output_file.sql
STEP 2 (in the same directory)
excute the following
psql -U postgresql_database_Username -d Database_name -W -f your_output_file.sql
your_output_file.sql is the file generated in step 1.
importing raster file to postgis is one time process n den u can use this raster table generated after importing ur tif file from the database for ur future raster queries
after this is done
u can use the postgis function (ST_Value) to get the pixel value of the raster u imported
*
SELECT ST_VALUE(e.rast, ST_SetSRID(St_MakePoint(Your_Longitude, Your_Latitude ), 4326))
FROM Your_Imported_raster_table_Name e
WHERE ST_Intersects(e.rast, ST_SetSRID(St_MakePoint(Your_Longitude, Your_Latitude), 4326));
*
Hope this will help you..
i would suggest that u go for 2nd method.
Best Answer
There are several ways to accomplish what you want. You've left out some specifics such as what extensions you have available and what license level you're using, as well as how thoroughly you want to sample the raster (ie how fine is your grid). I'll make some assumptions.
The simplest, if you want a point for every cell in the DEM, is the Raster to Point tool. A point will be created at the center of every cell with an attribute that is the value of that cell. There's also the Raster to ASCII or Raster to Float tools if those would produce formats more readily imported into ACAD than a shapefile (ACAD links in the last paragraph would indicate they would work just as well as the original DEM).
If you don't want every cell, you first need to use Create Fishnet to generate the grid. There's an option called Create Label Points, which will create a separate point feature class with points at the center of every grid cell. It's simplest to use these rather than all the extra steps needed to convert your grid lines to points at the intersections (which can be done, and is discussed in other questions). Since the points are at the center of the cells, you'll need to shift your grid origin -x and -y one half of the cell x/y value to get the points where the intersections would be. For example, 10m square grid, origin at -5,-5 will put the label point at 0,0. The grid can be any dimensions you like, including identical to the DEM cell size though that would just duplicate the first method. Once you have the points, you can use the Extract Values to Points tool from Spatial Analyst as you mention to get the elevation values. You could also use the Sample tool to generate a table (not shapefile) with x and y coordinate columns as well as an elevation column.
You can also convert the raster directly to a TIN using Raster to TIN if you have 3D Analyst, however I don't know how compatible Arc's TIN format is with ACAD (if at all). That said, as Aaron mentions you should be able to add a raster DEM directly to ACAD. I found this ACAD KB article that suggests you can. There's also this blog post and a helpfile / user's guide page that covers it. In fact a quick Google search turned up a lot of results on the subject.