Changes between Version 1 and Version 2 of DevWikiSpatialCollectionTutorial


Ignore:
Timestamp:
Sep 16, 2011, 8:51:38 PM (10 years ago)
Author:
bnordgren
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • DevWikiSpatialCollectionTutorial

    v1 v2  
    7676The salient point here is that value logic is separated from spatial logic. Different methods of determining values may be implemented and passed to the relation operator. This allows the relation operator code to be reused in other contexts (perhaps as a different function exposed to the user, perhaps as a utility object performing some supporting computation.) Likewise, the first-value evaluator may be reused as a component of any spatial collection where that behavior is appropriate.
    7777
     78In addition, `relation_op` is defined entirely in terms of the spatial collection interface. Therefore, the exact same code may be reused for other combinations of spatial data types: vector/raster; raster/vector; or vector/vector. The only difference is the very first step when the source data is initially wrapped into the spatial collection interface. Instead of four variants on roughly the same theme, there is ''one'' implementation to write and maintain.
     79
    7880Each of these spatial collections carries with it a `GBOX` named `extent`. When `relation_op` was instantiated, an approximate (or expected) extent was calculated based on the specified relation (intersection, union...) and the extents of the two inputs. It is important to point out that at this point, we have ''set up'' an operator, but we have not performed any actual computations. The `extent` of `relation_op` may be too big. It is best to take this with a grain of salt.
    7981
     82
     83== Creating an empty result raster ==
     84
     85When 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
     115Here 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
     117The 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
     119Finally, 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
     123Now 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
     125Again, 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
     127In 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
     163This 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
     165The 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
     169This 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
     178We 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.