Python Ctypes Examples
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()
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.
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)