... to extend a little bit presented GDAL functionality in my previous post I’ve decided to show you how to create raster file from list of coordinates stored in ASCII file. The intention is to create raster with pixel value taken from third coordinate. This can be anything from height, temperature, moisture, wind speed,...
I will show you how to to perform that with three tools I’m currently using GDAL, FME and ArcGIS. Unfortunately, this time I will limit myself only to these three tools but maybe in some post in the future I’ll include PostGIS, Oracle...
This tasks for GIS professionals is not demanding but many people with limited GIS knowledge need this and that is why I have decided to write down something about it or better to say demonstrate how to perform this task using three well known GIS tools.
The source ASCII file contains three columns with coordinates: easting northing and height.
This is segment of sample dataset (dtm_predicted_tides.csv):
445058.00, 6707614.00, 114.20 445060.00, 6707614.00, 114.22 445062.00, 6707614.00, 114.24 445064.00, 6707614.00, 114.25 445022.00, 6707616.00, 113.94 445024.00, 6707616.00, 113.96
... lets start with GDAL.
In my previous post I’ve used GDAL to generate concave hull.
This procedure very similar. Firstly we need to create vrt file (dem.vrt).
<OGRVRTDataSource>
<OGRVRTLayer name="test_dataset">
<LayerSRS>EPSG:23031</LayerSRS>
<SrcDataSource>dtm_predicted_tides.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
</OGRVRTLayer>
</OGRVRTDataSource>
To get extent use ogrinfo command. Parameters -al and -so will give you summary information for all layers (we have only one layer defined in vrt).
M:\__gdal_test>ogrinfo -al -so dem.vrt
INFO: Open of `dem.vrt'
using driver `VRT' successful.
Layer name: dtm_predicted_tides
Geometry: Point
Feature Count: 433752
Extent: (444448.000000, 6707606.000000) - (447314.000000, 6709738.000000)
Layer SRS WKT:
PROJCS["ED50 / UTM zone 31N",
GEOGCS["ED50",
DATUM["European_Datum_1950",
SPHEROID["International 1924",6378388,297,
AUTHORITY["EPSG","7022"]],
TOWGS84[-87,-98,-121,0,0,0,0],
AUTHORITY["EPSG","6230"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4230"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",3],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","23031"]]
field_1: String (0.0)
field_2: String (0.0)
field_3: String (0.0)
From data extent and resolution of data stored in ASCII file calculate dimension and extent for raster file you are about to generate. To get points centered in relation to pixel you should expand raster area for half resolution value. In this case resolution is the same for easting and northing axes. In this simple calculation do not forget that number of pixels (ie. raster resolution) can only be integer value.
rmaxx=maxx+resx/2 rmaxy=maxy+resx/2 rminx=minx-resy/2 rminy=miny-resy/2 dimx=(rmaxx-rminx)/resx dimy=(rmaxy-rminy)/resy txe: rminx rmaxx tye: rminy rmaxy outsize: dimx dimyCommand gdal_grid will produce desired grid file (dem.tiff). Search ellipse radius 1 and 2 should be smaller than data resolution.
gdal_grid -a nearest:radius1=1:radius2=1:nodata:-999
-txe 444447 447315 -tye 6707605 6709739
-outsize 1434 1067
-of GTiff -ot Float32
-l dtm_predicted_tides dem.vrt dtm_predicted_tides dem.tiff
To check results you can use gdalinfo command.
M:\__gdal_test>gdalinfo dem.tiff
Driver: GTiff/GeoTIFF
Files: dem.tiff
Size is 1434, 1067
Coordinate System is:
PROJCS["ED50 / UTM zone 31N",
GEOGCS["ED50",
DATUM["European_Datum_1950",
SPHEROID["International 1924",6378388,297.0000000000014,
AUTHORITY["EPSG","7022"]],
TOWGS84[-87,-98,-121,0,0,0,0],
AUTHORITY["EPSG","6230"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4230"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",3],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","23031"]]
Origin = (444447.000000000000000,6707605.000000000000000)
Pixel Size = (2.000000000000000,2.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 444447.000, 6707605.000) ( 1d59'19.59"E, 60d29'57.52"N)
Lower Left ( 444447.000, 6709739.000) ( 1d59'17.44"E, 60d31' 6.48"N)
Upper Right ( 447315.000, 6707605.000) ( 2d 2'27.50"E, 60d29'58.91"N)
Lower Right ( 447315.000, 6709739.000) ( 2d 2'25.46"E, 60d31' 7.87"N)
Center ( 445881.000, 6708672.000) ( 2d 0'52.50"E, 60d30'32.70"N)
Band 1 Block=1434x1 Type=Float32, ColorInterp=Gray
... here is FME example ...
This is even more simple. The first thing you have to define is csv data reader. In the way this is just like vrt file. For large datasets do not forget to limit number of scanned lines.
After that create geometry features form attributes using 3DPointReplacer
Use NumericRasterizer to create raster. You can define size of raster by entering number of cells by entering resolution (size of pixel) just like shown on the next figure.
In this moment everything is ready for to be written to some raster/grid format. In this example I’ve used Geotiff. On the writer you can add coordinate system definition.
Here is workflow for the whole process
...here is ArcGIS example...
In the Catalog Tree navigate to the 3D Analyst toolbox by expanding Toolboxes / System Toolboxes / 3D Analyst Tools. Double click ASCII 3D to Feature Class tool.
Enter all necessary values. That will create point features from the ASCII file. After that just start second command (see next figure) and the task is finished... for large dataset consider to drink a coffee or take a power nap :)
As I consider myself FME addict :) I just love model builder available in ArcGIS. Let me demonstrate how to create previous process using Model builder. Just open ModelBuilder window from Geoprocessing dropdown menu, drag tool from toolbox to Model builder and configure it just like it was demonstrated in previous steps. Process can be started from Model menu. Here is final look of process and created raster...
...and one more interesting functionality from ESRI... you can export this process to python script and then you have basic process coded in python and you can modify it and make it more powerful. It is very good start if you want to automatise some processes. Here is generated python code.
# ---------------------------------------------------------------------------
# process.py
# Created on: 2012-09-19 13:25:05.00000
# (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
# Check out any necessary licenses
arcpy.CheckOutExtension("3D")
# Local variables:
dtm_predicted_tides_csv__2_ = "M:\\__gdal_test\\dtm_predicted_tides.csv"
dtm_predicted_tides_MB = "\\\\...\\...\\gis_data.gdb\\dtm_predicted_tides_MB"
dtm_predicted_tides_MBr = "\\\\...\\...\\gis_data.gdb\\dtm_predicted_tides_MBr"
# Process: ASCII 3D to Feature Class
arcpy.ASCII3DToFeatureClass_3d("M:\\__gdal_test\\dtm_predicted_tides.csv", "XYZ", dtm_predicted_tides_MB, "POINT", "", "PROJCS['ED_1950_UTM_Zone_31N',GEOGCS['GCS_European_1950',DATUM['D_European_1950',SPHEROID['International_1924',6378388.0,297.0]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',3.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]", "", "", "DECIMAL_POINT")
# Process: Point to Raster
arcpy.PointToRaster_conversion(dtm_predicted_tides_MB, "Shape.Z", dtm_predicted_tides_MBr, "MAXIMUM", "NONE", "2")
Now just imagine how simple it will be to automatise this process for multiple files stored in some folder on filesystem.... with the help of little bit of python coding run trough folders, get file names, create list of filenames, create loop using generated python code and ... voilà ... you have new automated process :-)
....... I just want to mention here that FME and GDAL are also good integrated with python.
I hope that this post helps some of GIS beginners or some people that use GIS only occasionally....
....enjoy!
Thanks! Good luck with the job :)
ReplyDelete