| | 82 | |
| | 83 | == Creating an empty result raster == |
| | 84 | |
| | 85 | When creating a "raster result", one is typically talking about evaluating an expression at each grid cell in the raster which is to be returned. In the previous section, the creation of the `relation_op` operator set up the desired expression, but we have not yet performed any computations. We now need to define a raster which will contain the result, and then ask `relation_op` for a value over each of the grid cells. |
| | 86 | |
| | 87 | {{{ |
| | 88 | #!C |
| | 89 | /* make an empty raster to store the result */ |
| | 90 | result = rt_raster_new(res_width, res_height) ; |
| | 91 | |
| | 92 | /* copy srid and geo transform (except for offsets) */ |
| | 93 | rt_raster_set_srid(result, r1_pg->srid) ; |
| | 94 | rt_raster_set_scales(result, r1_pg->scaleX, r1_pg->scaleY) ; |
| | 95 | rt_raster_set_skews(result, r1_pg->skewX, r1_pg->skewY) ; |
| | 96 | |
| | 97 | /* compute new offsets because raster may be rotated |
| | 98 | * w.r.t. the extent. |
| | 99 | */ |
| | 100 | fit_raster_to_extent(res_extent, result) ; |
| | 101 | |
| | 102 | /* finally, wrap the relation operator to allow access by |
| | 103 | * indices of the result raster. |
| | 104 | */ |
| | 105 | relation_op_aligned = |
| | 106 | sc_create_raster_aligned_collection(relation_op, result) ; |
| | 107 | |
| | 108 | if (relation_op_aligned != NULL) { |
| | 109 | /* sample the relation operator into the raster */ |
| | 110 | sc_sampling_engine(relation_op_aligned, result, NULL) ; |
| | 111 | } |
| | 112 | |
| | 113 | }}} |
| | 114 | |
| | 115 | Here we create a new raster with a width and height computed to maintain the same x and y resolution as r1_pg, fitting within the expected extent of the result (`relation_op->extent`). We then copy all of the geolocation information except for the offsets into the result from `pg_r1`. This means that the pixels will have the same x and y scale and the same rotation as `pg_r1`. There is no particular reason that this must be the case: `pg_r1` is merely a convenient source of resolution and rotation information which is likely to be of interest to the user. Any geotransform parameters will do, as long as `result` is big enough to cover the expected result area and is in the same projection/srid. The function could even be defined to accept an empty raster from the user, if they want control over the output. |
| | 116 | |
| | 117 | The last step, `fit_raster_to_extent`, calculates and sets appropriate values for x and y offsets which position `result` in the extent of `relation_op`. With the location and alignment of `result` finalized, we then make it possible to query `relation_op` by the indices (pixel locations) of `result`. The constructor `sc_create_raster_aligned_collection` will add the `evaluateIndex` and `includesIndex` methods to any spatial collection. Again, this could be any spatial collection, including one which is derived from a purely geometric data type. The affect of this is to align the spatial collection to the grid defined by the provided raster. |
| | 118 | |
| | 119 | Finally, we get to the meat of any and all operations which return a raster filled with values: generating a value for every pixel with `sc_sampling_engine`. This is the topic of the next section. Note that in this section, we defined the geometry, but the raster itself is empty (meaning that it has no bands defined.) |
| | 120 | |
| | 121 | == Sampling a spatial collection == |
| | 122 | |
| | 123 | Now we have a way to generate values and we have a grid of points. These two come together in `sc_sampling_engine`. This function will create the appropriate number of new bands if necessary, but it can use existing bands if they are present. The only constraint is that the number of bands in the result raster must be the same as the number of elements returned by the spatial collection when it is asked for a value. |
| | 124 | |
| | 125 | Again, this code is written entirely against the spatial collection interface. This same sampling engine could be re-used to "rasterize" any spatial collection, including a wrapper around a single geometry which only returns two values: inside or outside. |
| | 126 | |
| | 127 | In the preliminaries, we have checked to ensure that the proper number of bands exist to receive data (or we have created them). We have also initialized a nodata ''vector'' (not a single nodata value). This is necessary because if the collection does not include a particular point, each band of the raster must still have a value set. The nodata vector allows for the possibility that different bands might have different nodata values. After these preliminaries, we have : |
| | 128 | |
| | 129 | {{{ |
| | 130 | #!C |
| | 131 | /* now sample the raster */ |
| | 132 | sample_pt = lwpoint_make2d(SRID_UNKNOWN, 0,0) ; |
| | 133 | for (j=0; j<height; j++) { |
| | 134 | for (i=0; i<width; i++) { |
| | 135 | int band_num ; |
| | 136 | VALUE *store_val ; |
| | 137 | |
| | 138 | /* set up the sample point */ |
| | 139 | sample_pt_p4d.x = i ; |
| | 140 | sample_pt_p4d.y = j ; |
| | 141 | ptarray_set_point4d(sample_pt->point,0,&sample_pt_p4d) ; |
| | 142 | |
| | 143 | /* evaluate the collection */ |
| | 144 | collection_val = sc_evaluateIndex(source, sample_pt) ; |
| | 145 | |
| | 146 | /* store the returned value or the nodata vector */ |
| | 147 | store_val = collection_val ; |
| | 148 | if (collection_val == NULL) { |
| | 149 | store_val = nodata_val ; |
| | 150 | } |
| | 151 | |
| | 152 | /* copy the values to the raster */ |
| | 153 | for (band_num=0; band_num < coll_bands; band_num++) { |
| | 154 | rt_band band ; |
| | 155 | |
| | 156 | band = rt_raster_get_band(result, band_num) ; |
| | 157 | rt_band_set_pixel(band, i, j, store_val->data[band_num]) ; |
| | 158 | } |
| | 159 | } |
| | 160 | } |
| | 161 | }}} |
| | 162 | |
| | 163 | This function, `sc_sampling_engine` may be reused whenever it is necessary to populate a raster with values from a spatial collection. It may be used with ''any'' spatial collection. In this case, the spatial collection is a "two input" spatial collection which performs an operation on the inputs. It could also be used on a spatial collection which represents only a single input. It is also possible to use `sc_sampling_engine` to populate a raster with the result of an equation (e.g., something which does not ultimately based on either a vector or raster data type. |
| | 164 | |
| | 165 | The salient point here is that `sc_sampling_engine` does not care how the values for each pixel are generated, nor does it care what those values are. It merely asks "What is the value at point '''P'''?" for every grid cell. It doesn't even care what algorithm was used to determine the size, shape, resolution or rotation of the raster. |
| | 166 | |
| | 167 | == Potential directions == |
| | 168 | |
| | 169 | This section lists some ideas about potentially useful spatial collection implementations: |
| | 170 | |
| | 171 | * Implement a spatial collection which evaluates an expression, plugging in values from a single wrapped collection. (one input mapalgebra) |
| | 172 | * Implement a spatial collection which evaluates an expression, plugging in values from two wrapped collections. (two input mapalgebra) |
| | 173 | * Currently, single geometries and single rasters are the only spatial objects which can be wrapped. Implement "overviews" which can select a set of geometries or rasters from a table based on some criteria from the other columns in the table. |
| | 174 | * The current implementation is very much geared towards evaluating data at a point, which is exactly what is needed to populate a raster with values. To generate "geometry" results, it will be necessary to add an iterator over all the geometric objects in the collection to the interface. Both vector and raster data can also answer the question "What spatial objects do you contain?" |
| | 175 | |
| | 176 | == Summary == |
| | 177 | |
| | 178 | We have gone a long way to erasing the distinction between raster and vector data by articulating two questions which both data types can answer. From these two questions, a set of interfaces have been designed. Implementing code against these interfaces instead of against the base data types allows the same algorithms to be applied to different data. This in turn promotes reuse and reduces the number of combinations which must be explicitly coded and maintained. |
