t.vect.observe.strds.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3. ############################################################################
  4. #
  5. # MODULE: t.vect.observe.strds
  6. # AUTHOR(S): Soeren Gebbert
  7. #
  8. # PURPOSE: Observe specific locations in a space time raster dataset over a period of time using vector points
  9. # COPYRIGHT: (C) 2011-2014 by the GRASS Development Team
  10. #
  11. # This program is free software under the GNU General Public
  12. # License (version 2). Read the file COPYING that comes with GRASS
  13. # for details.
  14. #
  15. #############################################################################
  16. #%module
  17. #% description: Observes specific locations in a space time raster dataset over a period of time using vector points.
  18. #% keyword: temporal
  19. #% keyword: sampling
  20. #% keyword: vector
  21. #% keyword: time
  22. #%end
  23. #%option G_OPT_V_INPUT
  24. #%end
  25. #%option G_OPT_STRDS_INPUTS
  26. #% key: strds
  27. #%end
  28. #%option G_OPT_STVDS_OUTPUT
  29. #%end
  30. #%option G_OPT_V_OUTPUT
  31. #% key: vector_output
  32. #% description: Name of the new created vector map that stores the sampled values in different layers
  33. #%end
  34. #%option
  35. #% key: columns
  36. #% type: string
  37. #% description: Names of the vector columns to be created and to store sampled raster values, one name for each STRDS
  38. #% required: yes
  39. #% multiple: yes
  40. #%end
  41. #%option G_OPT_DB_WHERE
  42. #%end
  43. #%option G_OPT_T_WHERE
  44. #% key: t_where
  45. #%end
  46. import grass.script as grass
  47. import grass.temporal as tgis
  48. import grass.script.raster as raster
  49. from grass.exceptions import CalledModuleError
  50. ############################################################################
  51. class Sample(object):
  52. def __init__(self, start = None, end = None, raster_names = None):
  53. self.start = start
  54. self.end = end
  55. if raster_names != None:
  56. self.raster_names = raster_names
  57. else:
  58. self.raster_names = []
  59. def __str__(self):
  60. return "Start: %s\nEnd: %s\nNames: %s\n"%(str(self.start), str(self.end), str(self.raster_names))
  61. ############################################################################
  62. def main():
  63. # Get the options
  64. input = options["input"]
  65. output = options["output"]
  66. vector_output = options["vector_output"]
  67. strds = options["strds"]
  68. where = options["where"]
  69. columns = options["columns"]
  70. tempwhere = options["t_where"]
  71. if where == "" or where == " " or where == "\n":
  72. where = None
  73. overwrite = grass.overwrite()
  74. # Check the number of sample strds and the number of columns
  75. strds_names = strds.split(",")
  76. column_names = columns.split(",")
  77. if len(strds_names) != len(column_names):
  78. grass.fatal(_("The number of columns must be equal to the number of space time raster datasets"))
  79. # Make sure the temporal database exists
  80. tgis.init()
  81. # We need a database interface
  82. dbif = tgis.SQLDatabaseInterfaceConnection()
  83. dbif.connect()
  84. mapset = grass.gisenv()["MAPSET"]
  85. out_sp = tgis.check_new_stds(output, "stvds", dbif, overwrite)
  86. samples = []
  87. first_strds = tgis.open_old_stds(strds_names[0], "strds", dbif)
  88. # Single space time raster dataset
  89. if len(strds_names) == 1:
  90. rows = first_strds.get_registered_maps(
  91. "name,mapset,start_time,end_time", tempwhere, "start_time", dbif)
  92. if not rows:
  93. dbif.close()
  94. grass.fatal(_("Space time raster dataset <%s> is empty") %
  95. out_sp.get_id())
  96. for row in rows:
  97. start = row["start_time"]
  98. end = row["end_time"]
  99. raster_maps = [row["name"] + "@" + row["mapset"],]
  100. s = Sample(start, end, raster_maps)
  101. samples.append(s)
  102. else:
  103. # Multiple space time raster datasets
  104. for name in strds_names[1:]:
  105. dataset = tgis.open_old_stds(name, "strds", dbif)
  106. if dataset.get_temporal_type() != first_strds.get_temporal_type():
  107. grass.fatal(_("Temporal type of space time raster datasets must be equal\n"
  108. "<%(a)s> of type %(type_a)s do not match <%(b)s> of type %(type_b)s"%\
  109. {"a":first_strds.get_id(),
  110. "type_a":first_strds.get_temporal_type(),
  111. "b":dataset.get_id(),
  112. "type_b":dataset.get_temporal_type()}))
  113. mapmatrizes = tgis.sample_stds_by_stds_topology("strds", "strds", strds_names,
  114. strds_names[0], False, None,
  115. "equal", False, False)
  116. for i in xrange(len(mapmatrizes[0])):
  117. isvalid = True
  118. mapname_list = []
  119. for mapmatrix in mapmatrizes:
  120. entry = mapmatrix[i]
  121. if entry["samples"]:
  122. sample = entry["samples"][0]
  123. name = sample.get_id()
  124. if name is None:
  125. isvalid = False
  126. break
  127. else:
  128. mapname_list.append(name)
  129. if isvalid:
  130. entry = mapmatrizes[0][i]
  131. map = entry["granule"]
  132. start, end = map.get_temporal_extent_as_tuple()
  133. s = Sample(start, end, mapname_list)
  134. samples.append(s)
  135. num_samples = len(samples)
  136. # Get the layer and database connections of the input vector
  137. vector_db = grass.vector.vector_db(input)
  138. # We copy the vector table and create the new layers
  139. if vector_db:
  140. # Use the first layer to copy the categories from
  141. layers = "1,"
  142. else:
  143. layers = ""
  144. first = True
  145. for layer in range(num_samples):
  146. layer += 1
  147. # Skip existing layer
  148. if vector_db and layer in vector_db and \
  149. vector_db[layer]["layer"] == layer:
  150. continue
  151. if first:
  152. layers += "%i" % (layer)
  153. first = False
  154. else:
  155. layers += ",%i" % (layer)
  156. vectmap = vector_output
  157. # We create a new vector map using the categories of the original map
  158. try:
  159. grass.run_command("v.category", input=input, layer=layers,
  160. output=vectmap, option="transfer",
  161. overwrite=overwrite)
  162. except CalledModuleError:
  163. grass.fatal(_("Unable to create new layers for vector map <%s>")
  164. % (vectmap))
  165. title = _("Observaion of space time raster dataset(s) <%s>") % (strds)
  166. description= _("Observation of space time raster dataset(s) <%s>"
  167. " with vector map <%s>") % (strds, input)
  168. # Create the output space time vector dataset
  169. out_sp = tgis.open_new_stds(output, "stvds",
  170. first_strds.get_temporal_type(),
  171. title, description,
  172. first_strds.get_semantic_type(),
  173. dbif, overwrite)
  174. dummy = out_sp.get_new_map_instance(None)
  175. # Sample the space time raster dataset with the vector
  176. # map at specific layer with v.what.rast
  177. count = 1
  178. for sample in samples:
  179. raster_names = sample.raster_names
  180. if len(raster_names) != len(column_names):
  181. grass.fatal(_("The number of raster maps in a granule must "
  182. "be equal to the number of column names"))
  183. # Create the columns creation string
  184. columns_string = ""
  185. for name, column in zip(raster_names, column_names):
  186. # The column is by default double precision
  187. coltype = "DOUBLE PRECISION"
  188. # Get raster map type
  189. raster_map = tgis.RasterDataset(name)
  190. raster_map.load()
  191. if raster_map.metadata.get_datatype() == "CELL":
  192. coltype = "INT"
  193. tmp_string = "%s %s,"%(column, coltype)
  194. columns_string += tmp_string
  195. # Remove last comma
  196. columns_string = columns_string[0:len(columns_string) - 1]
  197. # Try to add a column
  198. if vector_db and count in vector_db and vector_db[count]["table"]:
  199. try:
  200. grass.run_command("v.db.addcolumn", map=vectmap,
  201. layer=count, column=columns_string,
  202. overwrite=overwrite)
  203. except CalledModuleError:
  204. dbif.close()
  205. grass.fatal(_("Unable to add column %s to vector map <%s> "
  206. "with layer %i") % (columns_string, vectmap, count))
  207. else:
  208. # Try to add a new table
  209. grass.message("Add table to layer %i" % (count))
  210. try:
  211. grass.run_command("v.db.addtable", map=vectmap, layer=count,
  212. columns=columns_string, overwrite=overwrite)
  213. except CalledModuleError:
  214. dbif.close()
  215. grass.fatal(_("Unable to add table to vector map "
  216. "<%s> with layer %i") % (vectmap, count))
  217. # Call v.what.rast for each raster map
  218. for name, column in zip(raster_names, column_names):
  219. try:
  220. grass.run_command("v.what.rast", map=vectmap,
  221. layer=count, raster=name,
  222. column=column, where=where)
  223. except CalledModuleError:
  224. dbif.close()
  225. grass.fatal(_("Unable to run v.what.rast for vector map <%s> "
  226. "with layer %i and raster map <%s>") % \
  227. (vectmap, count, str(raster_names)))
  228. vect = out_sp.get_new_map_instance(dummy.build_id(vectmap,
  229. mapset, str(count)))
  230. vect.load()
  231. start = sample.start
  232. end = sample.end
  233. if out_sp.is_time_absolute():
  234. vect.set_absolute_time(start, end)
  235. else:
  236. vect.set_relative_time(
  237. start, end, first_strds.get_relative_time_unit())
  238. if vect.is_in_db(dbif):
  239. vect.update_all(dbif)
  240. else:
  241. vect.insert(dbif)
  242. out_sp.register_map(vect, dbif)
  243. count += 1
  244. out_sp.update_from_registered_maps(dbif)
  245. dbif.close()
  246. if __name__ == "__main__":
  247. options, flags = grass.parser()
  248. main()