source.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541
  1. import json
  2. import os
  3. import sys
  4. import uuid
  5. from ctypes import (
  6. addressof,
  7. byref,
  8. c_buffer,
  9. c_char_p,
  10. c_double,
  11. c_int,
  12. c_void_p,
  13. string_at,
  14. )
  15. from pathlib import Path
  16. from django.contrib.gis.gdal.driver import Driver
  17. from django.contrib.gis.gdal.error import GDALException
  18. from django.contrib.gis.gdal.prototypes import raster as capi
  19. from django.contrib.gis.gdal.raster.band import BandList
  20. from django.contrib.gis.gdal.raster.base import GDALRasterBase
  21. from django.contrib.gis.gdal.raster.const import (
  22. GDAL_RESAMPLE_ALGORITHMS,
  23. VSI_DELETE_BUFFER_ON_READ,
  24. VSI_FILESYSTEM_PREFIX,
  25. VSI_MEM_FILESYSTEM_BASE_PATH,
  26. VSI_TAKE_BUFFER_OWNERSHIP,
  27. )
  28. from django.contrib.gis.gdal.srs import SpatialReference, SRSException
  29. from django.contrib.gis.geometry import json_regex
  30. from django.utils.encoding import force_bytes, force_str
  31. from django.utils.functional import cached_property
  32. class TransformPoint(list):
  33. indices = {
  34. "origin": (0, 3),
  35. "scale": (1, 5),
  36. "skew": (2, 4),
  37. }
  38. def __init__(self, raster, prop):
  39. x = raster.geotransform[self.indices[prop][0]]
  40. y = raster.geotransform[self.indices[prop][1]]
  41. super().__init__([x, y])
  42. self._raster = raster
  43. self._prop = prop
  44. @property
  45. def x(self):
  46. return self[0]
  47. @x.setter
  48. def x(self, value):
  49. gtf = self._raster.geotransform
  50. gtf[self.indices[self._prop][0]] = value
  51. self._raster.geotransform = gtf
  52. @property
  53. def y(self):
  54. return self[1]
  55. @y.setter
  56. def y(self, value):
  57. gtf = self._raster.geotransform
  58. gtf[self.indices[self._prop][1]] = value
  59. self._raster.geotransform = gtf
  60. class GDALRaster(GDALRasterBase):
  61. """
  62. Wrap a raster GDAL Data Source object.
  63. """
  64. destructor = capi.close_ds
  65. def __init__(self, ds_input, write=False):
  66. self._write = 1 if write else 0
  67. Driver.ensure_registered()
  68. # Preprocess json inputs. This converts json strings to dictionaries,
  69. # which are parsed below the same way as direct dictionary inputs.
  70. if isinstance(ds_input, str) and json_regex.match(ds_input):
  71. ds_input = json.loads(ds_input)
  72. # If input is a valid file path, try setting file as source.
  73. if isinstance(ds_input, (str, Path)):
  74. ds_input = str(ds_input)
  75. if not ds_input.startswith(VSI_FILESYSTEM_PREFIX) and not os.path.exists(
  76. ds_input
  77. ):
  78. raise GDALException(
  79. 'Unable to read raster source input "%s".' % ds_input
  80. )
  81. try:
  82. # GDALOpen will auto-detect the data source type.
  83. self._ptr = capi.open_ds(force_bytes(ds_input), self._write)
  84. except GDALException as err:
  85. raise GDALException(
  86. 'Could not open the datasource at "{}" ({}).'.format(ds_input, err)
  87. )
  88. elif isinstance(ds_input, bytes):
  89. # Create a new raster in write mode.
  90. self._write = 1
  91. # Get size of buffer.
  92. size = sys.getsizeof(ds_input)
  93. # Pass data to ctypes, keeping a reference to the ctypes object so
  94. # that the vsimem file remains available until the GDALRaster is
  95. # deleted.
  96. self._ds_input = c_buffer(ds_input)
  97. # Create random name to reference in vsimem filesystem.
  98. vsi_path = os.path.join(VSI_MEM_FILESYSTEM_BASE_PATH, str(uuid.uuid4()))
  99. # Create vsimem file from buffer.
  100. capi.create_vsi_file_from_mem_buffer(
  101. force_bytes(vsi_path),
  102. byref(self._ds_input),
  103. size,
  104. VSI_TAKE_BUFFER_OWNERSHIP,
  105. )
  106. # Open the new vsimem file as a GDALRaster.
  107. try:
  108. self._ptr = capi.open_ds(force_bytes(vsi_path), self._write)
  109. except GDALException:
  110. # Remove the broken file from the VSI filesystem.
  111. capi.unlink_vsi_file(force_bytes(vsi_path))
  112. raise GDALException("Failed creating VSI raster from the input buffer.")
  113. elif isinstance(ds_input, dict):
  114. # A new raster needs to be created in write mode
  115. self._write = 1
  116. # Create driver (in memory by default)
  117. driver = Driver(ds_input.get("driver", "MEM"))
  118. # For out of memory drivers, check filename argument
  119. if driver.name != "MEM" and "name" not in ds_input:
  120. raise GDALException(
  121. 'Specify name for creation of raster with driver "{}".'.format(
  122. driver.name
  123. )
  124. )
  125. # Check if width and height where specified
  126. if "width" not in ds_input or "height" not in ds_input:
  127. raise GDALException(
  128. "Specify width and height attributes for JSON or dict input."
  129. )
  130. # Check if srid was specified
  131. if "srid" not in ds_input:
  132. raise GDALException("Specify srid for JSON or dict input.")
  133. # Create null terminated gdal options array.
  134. papsz_options = []
  135. for key, val in ds_input.get("papsz_options", {}).items():
  136. option = "{}={}".format(key, val)
  137. papsz_options.append(option.upper().encode())
  138. papsz_options.append(None)
  139. # Convert papszlist to ctypes array.
  140. papsz_options = (c_char_p * len(papsz_options))(*papsz_options)
  141. # Create GDAL Raster
  142. self._ptr = capi.create_ds(
  143. driver._ptr,
  144. force_bytes(ds_input.get("name", "")),
  145. ds_input["width"],
  146. ds_input["height"],
  147. ds_input.get("nr_of_bands", len(ds_input.get("bands", []))),
  148. ds_input.get("datatype", 6),
  149. byref(papsz_options),
  150. )
  151. # Set band data if provided
  152. for i, band_input in enumerate(ds_input.get("bands", [])):
  153. band = self.bands[i]
  154. if "nodata_value" in band_input:
  155. band.nodata_value = band_input["nodata_value"]
  156. # Instantiate band filled with nodata values if only
  157. # partial input data has been provided.
  158. if band.nodata_value is not None and (
  159. "data" not in band_input
  160. or "size" in band_input
  161. or "shape" in band_input
  162. ):
  163. band.data(data=(band.nodata_value,), shape=(1, 1))
  164. # Set band data values from input.
  165. band.data(
  166. data=band_input.get("data"),
  167. size=band_input.get("size"),
  168. shape=band_input.get("shape"),
  169. offset=band_input.get("offset"),
  170. )
  171. # Set SRID
  172. self.srs = ds_input.get("srid")
  173. # Set additional properties if provided
  174. if "origin" in ds_input:
  175. self.origin.x, self.origin.y = ds_input["origin"]
  176. if "scale" in ds_input:
  177. self.scale.x, self.scale.y = ds_input["scale"]
  178. if "skew" in ds_input:
  179. self.skew.x, self.skew.y = ds_input["skew"]
  180. elif isinstance(ds_input, c_void_p):
  181. # Instantiate the object using an existing pointer to a gdal raster.
  182. self._ptr = ds_input
  183. else:
  184. raise GDALException(
  185. 'Invalid data source input type: "{}".'.format(type(ds_input))
  186. )
  187. def __del__(self):
  188. if self.is_vsi_based:
  189. # Remove the temporary file from the VSI in-memory filesystem.
  190. capi.unlink_vsi_file(force_bytes(self.name))
  191. super().__del__()
  192. def __str__(self):
  193. return self.name
  194. def __repr__(self):
  195. """
  196. Short-hand representation because WKB may be very large.
  197. """
  198. return "<Raster object at %s>" % hex(addressof(self._ptr))
  199. def _flush(self):
  200. """
  201. Flush all data from memory into the source file if it exists.
  202. The data that needs flushing are geotransforms, coordinate systems,
  203. nodata_values and pixel values. This function will be called
  204. automatically wherever it is needed.
  205. """
  206. # Raise an Exception if the value is being changed in read mode.
  207. if not self._write:
  208. raise GDALException(
  209. "Raster needs to be opened in write mode to change values."
  210. )
  211. capi.flush_ds(self._ptr)
  212. @property
  213. def vsi_buffer(self):
  214. if not (
  215. self.is_vsi_based and self.name.startswith(VSI_MEM_FILESYSTEM_BASE_PATH)
  216. ):
  217. return None
  218. # Prepare an integer that will contain the buffer length.
  219. out_length = c_int()
  220. # Get the data using the vsi file name.
  221. dat = capi.get_mem_buffer_from_vsi_file(
  222. force_bytes(self.name),
  223. byref(out_length),
  224. VSI_DELETE_BUFFER_ON_READ,
  225. )
  226. # Read the full buffer pointer.
  227. return string_at(dat, out_length.value)
  228. @cached_property
  229. def is_vsi_based(self):
  230. return self._ptr and self.name.startswith(VSI_FILESYSTEM_PREFIX)
  231. @property
  232. def name(self):
  233. """
  234. Return the name of this raster. Corresponds to filename
  235. for file-based rasters.
  236. """
  237. return force_str(capi.get_ds_description(self._ptr))
  238. @cached_property
  239. def driver(self):
  240. """
  241. Return the GDAL Driver used for this raster.
  242. """
  243. ds_driver = capi.get_ds_driver(self._ptr)
  244. return Driver(ds_driver)
  245. @property
  246. def width(self):
  247. """
  248. Width (X axis) in pixels.
  249. """
  250. return capi.get_ds_xsize(self._ptr)
  251. @property
  252. def height(self):
  253. """
  254. Height (Y axis) in pixels.
  255. """
  256. return capi.get_ds_ysize(self._ptr)
  257. @property
  258. def srs(self):
  259. """
  260. Return the SpatialReference used in this GDALRaster.
  261. """
  262. try:
  263. wkt = capi.get_ds_projection_ref(self._ptr)
  264. if not wkt:
  265. return None
  266. return SpatialReference(wkt, srs_type="wkt")
  267. except SRSException:
  268. return None
  269. @srs.setter
  270. def srs(self, value):
  271. """
  272. Set the spatial reference used in this GDALRaster. The input can be
  273. a SpatialReference or any parameter accepted by the SpatialReference
  274. constructor.
  275. """
  276. if isinstance(value, SpatialReference):
  277. srs = value
  278. elif isinstance(value, (int, str)):
  279. srs = SpatialReference(value)
  280. else:
  281. raise ValueError("Could not create a SpatialReference from input.")
  282. capi.set_ds_projection_ref(self._ptr, srs.wkt.encode())
  283. self._flush()
  284. @property
  285. def srid(self):
  286. """
  287. Shortcut to access the srid of this GDALRaster.
  288. """
  289. return self.srs.srid
  290. @srid.setter
  291. def srid(self, value):
  292. """
  293. Shortcut to set this GDALRaster's srs from an srid.
  294. """
  295. self.srs = value
  296. @property
  297. def geotransform(self):
  298. """
  299. Return the geotransform of the data source.
  300. Return the default geotransform if it does not exist or has not been
  301. set previously. The default is [0.0, 1.0, 0.0, 0.0, 0.0, -1.0].
  302. """
  303. # Create empty ctypes double array for data
  304. gtf = (c_double * 6)()
  305. capi.get_ds_geotransform(self._ptr, byref(gtf))
  306. return list(gtf)
  307. @geotransform.setter
  308. def geotransform(self, values):
  309. "Set the geotransform for the data source."
  310. if len(values) != 6 or not all(isinstance(x, (int, float)) for x in values):
  311. raise ValueError("Geotransform must consist of 6 numeric values.")
  312. # Create ctypes double array with input and write data
  313. values = (c_double * 6)(*values)
  314. capi.set_ds_geotransform(self._ptr, byref(values))
  315. self._flush()
  316. @property
  317. def origin(self):
  318. """
  319. Coordinates of the raster origin.
  320. """
  321. return TransformPoint(self, "origin")
  322. @property
  323. def scale(self):
  324. """
  325. Pixel scale in units of the raster projection.
  326. """
  327. return TransformPoint(self, "scale")
  328. @property
  329. def skew(self):
  330. """
  331. Skew of pixels (rotation parameters).
  332. """
  333. return TransformPoint(self, "skew")
  334. @property
  335. def extent(self):
  336. """
  337. Return the extent as a 4-tuple (xmin, ymin, xmax, ymax).
  338. """
  339. # Calculate boundary values based on scale and size
  340. xval = self.origin.x + self.scale.x * self.width
  341. yval = self.origin.y + self.scale.y * self.height
  342. # Calculate min and max values
  343. xmin = min(xval, self.origin.x)
  344. xmax = max(xval, self.origin.x)
  345. ymin = min(yval, self.origin.y)
  346. ymax = max(yval, self.origin.y)
  347. return xmin, ymin, xmax, ymax
  348. @property
  349. def bands(self):
  350. return BandList(self)
  351. def warp(self, ds_input, resampling="NearestNeighbour", max_error=0.0):
  352. """
  353. Return a warped GDALRaster with the given input characteristics.
  354. The input is expected to be a dictionary containing the parameters
  355. of the target raster. Allowed values are width, height, SRID, origin,
  356. scale, skew, datatype, driver, and name (filename).
  357. By default, the warp functions keeps all parameters equal to the values
  358. of the original source raster. For the name of the target raster, the
  359. name of the source raster will be used and appended with
  360. _copy. + source_driver_name.
  361. In addition, the resampling algorithm can be specified with the "resampling"
  362. input parameter. The default is NearestNeighbor. For a list of all options
  363. consult the GDAL_RESAMPLE_ALGORITHMS constant.
  364. """
  365. # Get the parameters defining the geotransform, srid, and size of the raster
  366. ds_input.setdefault("width", self.width)
  367. ds_input.setdefault("height", self.height)
  368. ds_input.setdefault("srid", self.srs.srid)
  369. ds_input.setdefault("origin", self.origin)
  370. ds_input.setdefault("scale", self.scale)
  371. ds_input.setdefault("skew", self.skew)
  372. # Get the driver, name, and datatype of the target raster
  373. ds_input.setdefault("driver", self.driver.name)
  374. if "name" not in ds_input:
  375. ds_input["name"] = self.name + "_copy." + self.driver.name
  376. if "datatype" not in ds_input:
  377. ds_input["datatype"] = self.bands[0].datatype()
  378. # Instantiate raster bands filled with nodata values.
  379. ds_input["bands"] = [{"nodata_value": bnd.nodata_value} for bnd in self.bands]
  380. # Create target raster
  381. target = GDALRaster(ds_input, write=True)
  382. # Select resampling algorithm
  383. algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling]
  384. # Reproject image
  385. capi.reproject_image(
  386. self._ptr,
  387. self.srs.wkt.encode(),
  388. target._ptr,
  389. target.srs.wkt.encode(),
  390. algorithm,
  391. 0.0,
  392. max_error,
  393. c_void_p(),
  394. c_void_p(),
  395. c_void_p(),
  396. )
  397. # Make sure all data is written to file
  398. target._flush()
  399. return target
  400. def clone(self, name=None):
  401. """Return a clone of this GDALRaster."""
  402. if name:
  403. clone_name = name
  404. elif self.driver.name != "MEM":
  405. clone_name = self.name + "_copy." + self.driver.name
  406. else:
  407. clone_name = os.path.join(VSI_MEM_FILESYSTEM_BASE_PATH, str(uuid.uuid4()))
  408. return GDALRaster(
  409. capi.copy_ds(
  410. self.driver._ptr,
  411. force_bytes(clone_name),
  412. self._ptr,
  413. c_int(),
  414. c_char_p(),
  415. c_void_p(),
  416. c_void_p(),
  417. ),
  418. write=self._write,
  419. )
  420. def transform(
  421. self, srs, driver=None, name=None, resampling="NearestNeighbour", max_error=0.0
  422. ):
  423. """
  424. Return a copy of this raster reprojected into the given spatial
  425. reference system.
  426. """
  427. # Convert the resampling algorithm name into an algorithm id
  428. algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling]
  429. if isinstance(srs, SpatialReference):
  430. target_srs = srs
  431. elif isinstance(srs, (int, str)):
  432. target_srs = SpatialReference(srs)
  433. else:
  434. raise TypeError(
  435. "Transform only accepts SpatialReference, string, and integer "
  436. "objects."
  437. )
  438. if target_srs.srid == self.srid and (not driver or driver == self.driver.name):
  439. return self.clone(name)
  440. # Create warped virtual dataset in the target reference system
  441. target = capi.auto_create_warped_vrt(
  442. self._ptr,
  443. self.srs.wkt.encode(),
  444. target_srs.wkt.encode(),
  445. algorithm,
  446. max_error,
  447. c_void_p(),
  448. )
  449. target = GDALRaster(target)
  450. # Construct the target warp dictionary from the virtual raster
  451. data = {
  452. "srid": target_srs.srid,
  453. "width": target.width,
  454. "height": target.height,
  455. "origin": [target.origin.x, target.origin.y],
  456. "scale": [target.scale.x, target.scale.y],
  457. "skew": [target.skew.x, target.skew.y],
  458. }
  459. # Set the driver and filepath if provided
  460. if driver:
  461. data["driver"] = driver
  462. if name:
  463. data["name"] = name
  464. # Warp the raster into new srid
  465. return self.warp(data, resampling=resampling, max_error=max_error)
  466. @property
  467. def info(self):
  468. """
  469. Return information about this raster in a string format equivalent
  470. to the output of the gdalinfo command line utility.
  471. """
  472. return capi.get_ds_info(self.ptr, None).decode()