Difference between revisions of "Python Ctypes Examples"

From GRASS-Wiki
Jump to navigation Jump to search
(note what's grass7)
 
(8 intermediate revisions by 4 users not shown)
Line 1: Line 1:
See [[GRASS and Python]] for more information.
* See [[GRASS and Python]] for more information.


== Accessing feature geometry ==
* See also [http://trac.osgeo.org/grass/browser/grass/trunk/doc/python/raster_example_ctypes.py raster], [http://trac.osgeo.org/grass/browser/grass/trunk/doc/python/vector_example_ctypes.py vector] sample scripts. (GRASS 7)
 
 
 
 
__TOC__
== Raster ==
 
=== Working with (geographic) coordinates for each raster cell and the computational region ===
 
The ctypes interfaces are thin wrappers around the corresponding C
functions. Any script which uses them should have essentially the same
structure as a module written in C. Structures can be allocated by
using the structure name as a function, and passed by reference using
[http://docs.python.org/library/ctypes.html ctypes]' byref() function.
 
The second argument to G_col_to_easting() and G_row_to_northing() is a
"const struct Cell_head *". The information can be obtained from
G_get_window(), e.g.:
 
<source lang="python">
      import sys
      import grass.lib.gis as gis
      from ctypes import *
 
      gis.G_gisinit(sys.argv[0])
 
      region = gis.Cell_head()
      gis.G_get_window(byref(region))
      coor_col = gis.G_col_to_easting((thecolumn + 0.5), byref(region))
      coor_row = gis.G_row_to_northing((therow + 0.5), byref(region))
</source>
 
=== Basic raster metadata access methods ===
 
This is a simple demonstration script to show the ctypes style access to some of the raster metadata.
 
(''GRASS 7'')
 
<source lang="python">
from ctypes import *
import grass.lib.gis as libgis
import grass.lib.raster as libraster
 
import grass.script as grass
import grass.temporal as tgis
 
def print_raster_info(name):
 
  grass.run_command("r.mapcalc", expr=name + " = 5.0", overwrite=True)
  grass.run_command("r.timestamp", map=name, date="13 Jan 2003 / 15 Jan 2003")
 
  # Create a region struct
  region = libgis.Cell_head()
 
  # Read the raster region
  libraster.Rast_get_cellhd("test", "", byref(region))
  print region.cols
  print region.rows
  print region.north
  print region.south
  print region.east
  print region.west
  print region.ns_res
  print region.ew_res
  print region.top
  print region.bottom
  print region.tb_res
  print region.depths
  print region.proj
 
  # Get Map type
  maptype = libraster.Rast_map_type(name, "")
 
  if maptype == libraster.DCELL_TYPE:
    print "DCELL"
  elif maptype == libraster.FCELL_TYPE:
    print "FCELL"
  elif maptype == libraster.CELL_TYPE:
    print "CELL"
  else:
    print "Unknown type"
 
  # Read range
  if libraster.Rast_map_is_fp(name, ""):
    print "Floating point map"
    range = libraster.FPRange()
    libraster.Rast_read_fp_range(name, "", byref(range))
    min = libgis.DCELL()
    max = libgis.DCELL()
    libraster.Rast_get_fp_range_min_max(byref(range), byref(min), byref(max))
  else:
    print "Cell map"
    range = libraster.Range()
    libraster.Rast_read_range(name, "", byref(range))
    min = libgis.CELL()
    max = libgis.CELL()
    libraster.Rast_get_fp_range_min_max(byref(range), byref(min), byref(max))
   
  minimum = min.value
  maximum = max.value
  print minimum
  print maximum
 
  # Read time stamp
  tstmp = libgis.TimeStamp()
  libgis.G_read_raster_timestamp(name, "", byref(tstmp))
 
  print tstmp.count
 
  dt1 = tstmp.dt[0]
  dt2 = tstmp.dt[1]
 
  if tstmp.count >= 1:
    print dt1.mode
    print dt1.year
    print dt1.month
    print dt1.day
    print dt1.hour
    print dt1.minute
    print dt1.second
  if tstmp.count > 1:
    print dt2.mode
    print dt2.year
    print dt2.month
    print dt2.day
    print dt2.hour
    print dt2.minute
    print dt2.second
 
 
if __name__ == "__main__":
    print_raster_info("test")
</source>
 
== Vector ==
 
=== Accessing feature geometry ===


This script switches X, Y coordinates and multiple them by -1.
This script switches X, Y coordinates and multiple them by -1.


<source lang="python">
<source lang="python">
#!/usr/bin/python
import sys
import sys
import grass.script as grass


from grass.lib.vector import *
from grass.lib.vector import *
Line 28: Line 165:
nlines = Vect_get_num_lines(map_info)
nlines = Vect_get_num_lines(map_info)


for line in range(1, nlines):
for line in range(1, nlines + 1):
     ltype = Vect_read_line(map_info, points, cats, line)
     ltype = Vect_read_line(map_info, points, cats, line)
     if line % 100 == 0:
     if line % 100 == 0:
Line 41: Line 178:
                       points, cats)
                       points, cats)
     line += 1
     line += 1
sys.stderr.write("\r")


Vect_destroy_line_struct(points)
Vect_destroy_line_struct(points)
Vect_destroy_cats_struct(cats)
Vect_destroy_cats_struct(cats)


Vect_build_partial(map_info, GV_BUILD_NONE)
Vect_build(map_info)
Vect_close(map_info)
Vect_close(map_info)
</source>
</source>
== See also ==
* http://docs.python.org/library/ctypes.html


{{Python}}
{{Python}}

Latest revision as of 10:02, 20 May 2012



Raster

Working with (geographic) coordinates for each raster cell and the computational region

The ctypes interfaces are thin wrappers around the corresponding C functions. Any script which uses them should have essentially the same structure as a module written in C. Structures can be allocated by using the structure name as a function, and passed by reference using ctypes' byref() function.

The second argument to G_col_to_easting() and G_row_to_northing() is a "const struct Cell_head *". The information can be obtained from G_get_window(), e.g.:

       import sys
       import grass.lib.gis as gis
       from ctypes import *

       gis.G_gisinit(sys.argv[0])

       region = gis.Cell_head()
       gis.G_get_window(byref(region))
       coor_col = gis.G_col_to_easting((thecolumn + 0.5), byref(region))
       coor_row = gis.G_row_to_northing((therow + 0.5), byref(region))

Basic raster metadata access methods

This is a simple demonstration script to show the ctypes style access to some of the raster metadata.

(GRASS 7)

from ctypes import *
import grass.lib.gis as libgis
import grass.lib.raster as libraster

import grass.script as grass
import grass.temporal as tgis

def print_raster_info(name):
  
  grass.run_command("r.mapcalc", expr=name + " = 5.0", overwrite=True)
  grass.run_command("r.timestamp", map=name, date="13 Jan 2003 / 15 Jan 2003")
  
  # Create a region struct
  region = libgis.Cell_head()

  # Read the raster region
  libraster.Rast_get_cellhd("test", "", byref(region))
  print region.cols
  print region.rows
  print region.north
  print region.south
  print region.east
  print region.west
  print region.ns_res
  print region.ew_res
  print region.top
  print region.bottom
  print region.tb_res
  print region.depths
  print region.proj
  
  # Get Map type
  maptype = libraster.Rast_map_type(name, "")
  
  if maptype == libraster.DCELL_TYPE:
    print "DCELL"
  elif maptype == libraster.FCELL_TYPE:
    print "FCELL"
  elif maptype == libraster.CELL_TYPE:
    print "CELL"
  else:
    print "Unknown type"
  
  # Read range
  if libraster.Rast_map_is_fp(name, ""):
    print "Floating point map"
    range = libraster.FPRange()
    libraster.Rast_read_fp_range(name, "", byref(range))
    min = libgis.DCELL()
    max = libgis.DCELL()
    libraster.Rast_get_fp_range_min_max(byref(range), byref(min), byref(max))
  else:
    print "Cell map"
    range = libraster.Range()
    libraster.Rast_read_range(name, "", byref(range))
    min = libgis.CELL()
    max = libgis.CELL()
    libraster.Rast_get_fp_range_min_max(byref(range), byref(min), byref(max))
    
  minimum = min.value
  maximum = max.value
  print minimum
  print maximum
  
  # Read time stamp
  tstmp = libgis.TimeStamp()
  libgis.G_read_raster_timestamp(name, "", byref(tstmp))
  
  print tstmp.count
  
  dt1 = tstmp.dt[0]
  dt2 = tstmp.dt[1]
  
  if tstmp.count >= 1:
    print dt1.mode
    print dt1.year
    print dt1.month
    print dt1.day
    print dt1.hour
    print dt1.minute
    print dt1.second
  if tstmp.count > 1:
    print dt2.mode
    print dt2.year
    print dt2.month
    print dt2.day
    print dt2.hour
    print dt2.minute
    print dt2.second
  

if __name__ == "__main__":
    print_raster_info("test")

Vector

Accessing feature geometry

This script switches X, Y coordinates and multiple them by -1.

#!/usr/bin/python

import sys

from grass.lib.vector import *

if len(sys.argv) < 2:
    sys.exit("Usage: %s: vector" % sys.argv[0])

name = sys.argv[1]

map_info = pointer(Map_info())
points = Vect_new_line_struct()
cats = Vect_new_cats_struct()
c_points = points.contents

level = Vect_open_update(map_info, name, "")
if level < 2:
    sys.exit("Topology not available")

nlines = Vect_get_num_lines(map_info)

for line in range(1, nlines + 1):
    ltype = Vect_read_line(map_info, points, cats, line)
    if line % 100 == 0:
        sys.stderr.write("\r%d" % line)
    
    for i in range(c_points.n_points):
        x = c_points.x[i]
        c_points.x[i] = -1 * c_points.y[i]
        c_points.y[i] = -1 * x
    
    Vect_rewrite_line(map_info, line, ltype,
                      points, cats)
    line += 1

sys.stderr.write("\r")

Vect_destroy_line_struct(points)
Vect_destroy_cats_struct(cats)

Vect_build_partial(map_info, GV_BUILD_NONE)
Vect_build(map_info)
Vect_close(map_info)

See also