Difference between revisions of "Python Ctypes Examples"

From GRASS-Wiki
Jump to: navigation, search
(+raster example from Glynn)
(Raster)
Line 28: Line 28:
 
</source>
 
</source>
  
 +
=== Basic raster metadata access methods ===
 +
 +
This is a simple demonstration script to show the ctypes style access to some of the raster metadata.
 +
 +
<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 ==
 
== Vector ==

Revision as of 08:26, 7 April 2012

See GRASS and Python for more information. See also raster, vector sample scripts.

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()
       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.

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