Python Ctypes Examples
Jump to navigation
Jump to search
- See GRASS and Python for more information.
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)