Friday, September 21, 2012

PyDev Installation

keywords: Eclipse, python, pydev, installation

I have just installed Eclipse Version: Juno Release (Build id: 20120614-1722) and I’ve wanted to install Python IDE PyDev. When I’ve started with the installation procedure I’ve run into problem.... here are the error log and the screen shot ....

HTTP Proxy Authentication Required: http://pydev.org/updates/content.xml
HTTP Proxy Authentication Required: http://pydev.org/updates/content.xml
Proxy Authentication Required




Luckily the solution was simple. It was only required to add some configuration options to eclipse.ini file. Instructions how to solve this issue I’ve found on internet were not so precise and possibly could make some confusion ... so I’ve decided to make one short post on this topic.

.... locate your eclipse.ini file and add at the end these lines:
-Dorg.eclipse.ecf.provider.filetransfer.excludeContributors=org.eclipse.ecf.provider.filetransfer.httpclient
-Dhttp.proxyPort=8080
-Dhttp.proxyHost=myproxy
-Dhttp.proxyUser=mydomain\myusername
-Dhttp.proxyPassword=mypassword
-Dhttp.nonProxyHosts=localhost|127.0.0.1
It is important to put these lines to last block after -vmargs. Here is the whole ini file (added lines 22-27).
-startup
plugins/org.eclipse.equinox.launcher_1.3.0.v20120522-1813.jar
--launcher.library
plugins/org.eclipse.equinox.launcher.win32.win32.x86_1.1.200.v20120522-1813
-product
org.eclipse.epp.package.jee.product
--launcher.defaultAction
openFile
--launcher.XXMaxPermSize
256M
-showsplash
org.eclipse.platform
--launcher.XXMaxPermSize
256m
--launcher.defaultAction
openFile
-vmargs
-Dosgi.requiredJavaVersion=1.5
-Dhelp.lucene.tokenizer=standard
-Xms40m
-Xmx512m
-Dorg.eclipse.ecf.provider.filetransfer.excludeContributors=org.eclipse.ecf.provider.filetransfer.httpclient
-Dhttp.proxyPort=8080
-Dhttp.proxyHost=myproxy
-Dhttp.proxyUser=mydomain\myusername
-Dhttp.proxyPassword=mypassword
-Dhttp.nonProxyHosts=localhost|127.0.0.1

...and don’t forget to restart Eclipse :-)


....here are more details.

... enjoy ...

Wednesday, September 19, 2012

Create grid/raster from ASCII file using GDAL, FME and ArcGIS

keywords: ASCII file, raster, grid, GDAL, OGR, FME, ArcGIS, python

... 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 dimy
Command 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!

Friday, September 14, 2012

How to generate Concave Hull from regularly spaced point data using GDAL

keywords: GDAL, OGR, vrt, gdal_grid, gdal_contour, concave hull, grid/raster processing

If you have to to process data delivered in plain ASCII files you have probably run into concave hull problem. Contrary to convex hull problem GIS tools usually don’t have some builtin function that can help you with this issue. There is a good demonstration of PostGIS ST_ConcaveHull function

This is just to illustrate multiple possible solution to concave hull problem.

I this post presented scenario is not generic and it works best for regularly spaced points or almost regularly spaced points. With some adjustments can be applied for irregularly scattered points. For now I’ll describe solution for regularly spaced points. Although this procedure is not general I hope it will eventually help somebody.

The idea is to create raster/grid file from point data and to identify areas (pixels) that doesn’t contain any point. The easiest way is to produce grid in the way that value of each pixel represents number of points. For regularly spaced data if everything is setup well we should get only values 0 and 1.

From the created grid just generate contour with value 1. It represents concave hull for the regularly spaced data. If you have almost regularly spaced data you can apply the same principle and just generate contour with value 0.1 (this also works for regularly spacing) . To illustrate difference here is the figure with three contours with values 0.1, 0.5 and 1.

I’ll will start with simple ASCII file with data organized into three columns: x, y, z. That is very common data format that many professional use to share point data. Here is just few lines from used file dtm_predicted_tides.csv (in this scenario file does not contain headers).
446034.00, 6707606.00, 113.92
445036.00, 6707606.00, 113.95
445038.00, 6707606.00, 113.98
445040.00, 6707606.00, 114.00
...
First we have to create VRT (virtual format) file because GDAL cannot directly read and “understand” ASCII file. VRT file details you can find here.

For this application (no headers in ASCII file) VRT file will be (dem.vrt):
<OGRVRTDataSource>
    <OGRVRTLayer name="dtm_predicted_tides">
       <SrcDataSource>dtm_predicted_tides.csv</SrcDataSource>
               <GeometryType>wkbPoint</GeometryType>
               <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
    </OGRVRTLayer>
</OGRVRTDataSource>
To test it use command ogrinfo:
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:
(unknown)
field_1: String (0.0)
field_2: String (0.0)
field_3: String (0.0)
From this info data you will need extent to calculate some parameters that you will need in next step. To assure that generated contour line is closed it is necessary to extend area for at least one pixel value. The size of pixel is determined by resolution of original data. For regularly spaced points this value is equal to distance of two neighboring points along the easting and northing axes. For almost regularly spaced data you have to estimate this value bearing in mind the proposed scenario for concave hull estimation.

Taking into account generated extent and resolution of 2m it can be easily calculated that the generated grid should be sized 1433 x 1066. To extend area and to assure that point is located in the middle of generated pixel it is required to extend area for one and half pixel and in this case that is 3m.

Now we have all values necessary to run command. Notice that the size of grid is also changed because original extent was changed.
gdal_grid -a count:radius1=1:radius2=1
          -txe 444445 447317   -tye 6707603 6709741
          -outsize 1436 1069  
          -of GTiff -ot Float32 
          -l dtm_predicted_tides dem.vrt  concavehull_ext2.tiff
Usage: gdal_grid [--help-general] [--formats]
    [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
         CInt16/CInt32/CFloat32/CFloat64}]
    [-of format] [-co "NAME=VALUE"]
    [-zfield field_name]
    [-a_srs srs_def] [-spat xmin ymin xmax ymax]
    [-clipsrc |WKT|datasource|spat_extent]
    [-clipsrcsql sql_statement] [-clipsrclayer layer]
    [-clipsrcwhere expression]
    [-l layername]* [-where expression] [-sql select_statement]
    [-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]
    [-a algorithm[:parameter1=value1]*]    [-q]
    <src_datasource> <dst_filename>

Available algorithms and parameters with their's defaults:
    Inverse distance to a power (default)
invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0
    Moving average
       average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0
    Nearest neighbor
       nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0
    Various data metrics
       :radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0
       possible metrics are:
           minimum
           maximum
           range
           count
           average_distance
           average_distance_pts
Notice that we have used count metrics with the boath radius of search ellipse set to 1. The intention was to get one point per pixel. With data that are not regularly spaced this is not possible to get so in that case you should just play a little bit firstly with pixel size and then by search ellipse radius. At the end just interpolate... lets say 0.1 contour and you can get satisfactory results. .... Use gdalinfo to check generated raster
M:\__gdal_test>gdalinfo concavehull_ext2.tiff
Driver: GTiff/GeoTIFF
Files: concavehull_ext2.tiff
      concavehull_ext2.tiff.aux.xml
Size is 1436, 1069
Coordinate System is:
PROJCS["ED_1950_UTM_Zone_31N",
    GEOGCS["GCS_European_1950",
       DATUM["European_Datum_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],
    AUTHORITY["EPSG","23031"]]
Origin = (444445.000000000000000,6707603.000000000000000)
Pixel Size = (2.000000000000000,2.000000000000000)
Image Structure Metadata:
 INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  444445.000, 6707603.000) (  1d59'19.46"E, 60d29'57.45"N)
Lower Left  (  444445.000, 6709741.000) (  1d59'17.31"E, 60d31' 6.55"N)
Upper Right (  447317.000, 6707603.000) (  2d 2'27.63"E, 60d29'58.84"N)
Lower Right (  447317.000, 6709741.000) (  2d 2'25.59"E, 60d31' 7.94"N)
Center      (  445881.000, 6708672.000) (  2d 0'52.50"E, 60d30'32.70"N)
Band 1 Block=1436x1 Type=Float32, ColorInterp=Gray
 Min=0.000 Max=1.000
 Minimum=0.000, Maximum=1.000, Mean=0.283, StdDev=0.450
 Metadata:
    STATISTICS_MAXIMUM=1
    STATISTICS_MEAN=0,28255913031469
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0,45024393416031
Pixel size and corner coordinates should be as it was previously calculated. You won’t get info about statistics and coordinate system but for now that is not important.

In this moment we are ready to create concave hull or to be more precise we are able to generate contours. The following command will generate shape file with 0.1 contour and attribute elev that holds contour value.
gdal_contour -a elev concavehull_ext2.tiff concavehull_ext2_z.shp -fl 0.1
Usage: gdal_contour [-b ] [-a ] [-3d] [-inodata]
                   [-snodata n] [-f ] [-i ]
                   [-off ] [-fl  ...]
                   [-nln ] [-q]
                    
Here is generated concave hull for my case. Gray pixels contain no points and for the other pixels count value is equal to 1 (that was intention from the beginning). If you need you can generate multiple contours at once. In that case you can use elev attribute to extract what you intend to use for concave hull..
gdal_contour -a elev concavehull_ext2.tiff contours.shp -fl 0.1 0.5 0.7
If you need to visualise points you can use ogr2ogr to create shape file from vrt.
ogr2ogr dem.shp dem.vrt
Here is segment with metric raster, points and concave hull (contour=1). It is good practice to use coordinate system definition from the beginning of process. For that you should add one subelement to the vrt OGRVRTLayer element.
<LayerSRS>EPSG:23031</LayerSRS>
The previous is sufficient to propagate coordinate system definition to all generated files.

The whole process is very simple and can be summarized into few lines:
1. create vrt file and check it (ogrinfo); optionaly crate shape file ogr2ogr
2. create metric raster gdal_grid (check it with gdalinfo)
3. generate concave hull (gdal_contour)

... just play with your data and enjoy! :-)