Subscribe to this thread
Home - General / All posts - Creating a mesh/grid using raster cell boundaries in M9
gjsa31 post(s)
#11-Jul-18 02:48

Have been searching for ways to create a regular grid as vector drawing, but nothing so far that specifically applies in M9.

The links I've found so far:

http://www.georeference.org/doc/view_grid.htm

http://www.georeference.org/Forum/t132034.8

The trace areas feature is fantastic, but sometimes I just want the pixels turned into individual polygons regardless of their raster value. The equivalent ESRI process would be Create Fishnet.

Dimitri


4,899 post(s)
#11-Jul-18 06:07

You could do that in SQL or a script today. Point and click tools for grids are coming.

Grids as virtual display overlays and automatically creating grids as vectors have been in the wishlist for 9 for a while. They're just not as high a priority as other things, at least not in terms of what users have asked for. [...oh, no! ... not that link again! :-) ]

In my personal work I rarely need to create a regular grid, but if I do I just use 8 to create the thing and then use it in 9.

I have to admit to being curious:

sometimes I just want the pixels turned into individual polygons regardless of their raster value.

That's a very small grid in most cases, if I understand your desire as being to create a grid so small each cell in the grid is one pixel size. Why do you want to do that? If it is for operations on rasters, perhaps there is an alternative workflow.

gjsa31 post(s)
#11-Jul-18 07:55

It is a strange use case, and probably not one that would be common to other users.

Basically I need to systematically edit a raster at the pixel level - the easiest way to do that is to do vector operations (Overlay, adjacent - for example) on a grid that contains the values of individual pixels as polygons or perhaps centroids (points).

The final outcome would be all pixels of the same value merged into areas, but only after this editing step is completed.

For the time being I can use ArcGIS or I could buy M8 by the sounds of it.

tjhb
8,093 post(s)
#11-Jul-18 08:14

The first part of this—creating the pixel-aligned grid—is actually pretty easy in Manifold 9 SQL. Easier and mire natural than in Manifold 8.

But how about the next step, the systematic editing? And other steps.

The more you can say about the workflow you have in mind, the better. E.g. what you would do in Arc, as explained to soneone who does not have Arc, and paying attention to the purpose at each step.

Whatever you need to do will be possible, either now or soon (now I think). But with overhead—which it would be best to try to judge in advance.

gjsa31 post(s)
#13-Jul-18 08:01

Is there an example of a pixel-aligned vector grid creation in Manifold 9 SQL you could post when you get a chance?

One follow-up question also:

I've decided that extracting raster values to said grid is far more efficient than creating a vector from the raster and trying to overlay that onto the grid. But, pretty sure M9 doesn't do raster value extractions to shapes yet? If I'm wrong, happy to be corrected.

SAGA 6 does raster value extractions fast - but the slow part is exporting/importing very large shapefiles, so the more I can do in M9 the better.

Dimitri


4,899 post(s)
#13-Jul-18 08:24

But, pretty sure M9 doesn't do raster value extractions to shapes yet?

Yes, it does.

See the Example: Trace Vector Areas from Raster Pixels topic for an example.

tjhb
8,093 post(s)
#13-Jul-18 10:05
Is there an example of a pixel-aligned vector grid creation in Manifold 9 SQL you could post when you get a chance?

No. You need to say what you need to do with the grid first, as I said.

One obvious reason for that: you have said you need a pixel-aligned grid, that is one thing. But I suspect you actually mean a grid of pixel cells along with pixel values, which is different.

So to avoid a wild goose chase...

gjsa31 post(s)
#16-Jul-18 14:35

A long series of SQL commands I am writing using M9 templates will follow soon when they are complete.

They all begin with some input layers: one of which is a vector grid of raster pixel values and also an empty vector grid with the same dimensions and alignment. If I could create the grid of pixel values direct from the raster the second empty grid is easy (just need to copy and nullify all the values).

I appreciate the need to fully understand what I am trying to achieve, and I will complete the picture soon.

For now just to reiterate that Trace Areas is great and I use it often, but in this case I do need the pixels vectorised individually and haven't figured out how in M9.

tjhb
8,093 post(s)
#16-Jul-18 20:39

That's definitely clear enough, thanks. I'll see if I can get it right. As you'll see it's not hard, lots of built-in tools to make it straightforward.

mdsumner


4,204 post(s)
#16-Jul-18 23:29

Interesting to me is the abstractability of this, the source is six numbers (offset, scale, dimension) and always the code to retrieve these seems harder work than the explicit polygon generation. Obviously there's a lot of pipeline potential, keeping things implicit, the data transfer is trivial with a single pixel index.


https://github.com/mdsumner

tjhb
8,093 post(s)
#17-Jul-18 07:17

Here's a first go.

A few things to note.

First, this requires Manifold 9 >= 9.0.167.6, released today.

Secondly, the source data type is hard coded (here INT32X3, e.g. standard RGB or RGB). This needs to be adjusted (in three places) for other data types. I think that's currently unavoidable (unless we wrap in a script).

Third, there is a Rect filter, commented out. So all pixels from all tiles are returned, including those falling outside the image Rect. Uncomment the filter to remove peripheral cells for pixels that don't show.

Lastly, I don't think there is anything here that would benefit from multiple threads, but I haven't tested that yet.

Execution time isn't bad. For an image 5761 x 8505px, I get a time of 16 minutes on this i7 laptop, producing 49,397,760 cells (pixel areas).

--SQL9

-- source/target vector type is hardcoded (see notes)

-- Rect filter inactive (see WHERE clause)

--

-- source table and image

VALUE @table TABLE = [BA32_GRIDLESS_GeoTifv1-06 Tiles];

VALUE @image TABLE = [BA32_GRIDLESS_GeoTifv1-06]-- NB image qua TABLE

--

-- image metadata

VALUE @field NVARCHAR = ComponentProperty(@image'FieldTile');

VALUE @dimXY INT32X2 = CAST(ComponentProperty(@table'FieldTileSize.' + @field) AS INT32X2);

VALUE @dimX INT32 = VectorValue(@dimXY, 0);

VALUE @dimY INT32 = VectorValue(@dimXY, 1);

VALUE @rect INT32X4 = CAST(ComponentProperty(@image'Rect'AS INT32X4);

VALUE @minX INT32 = CAST(VectorValue(@rect, 0) AS INT32);

VALUE @minY INT32 = CAST(VectorValue(@rect, 1) AS INT32);

VALUE @maxX INT32 = CAST(VectorValue(@rect, 2) AS INT32);

VALUE @maxY INT32 = CAST(VectorValue(@rect, 3) AS INT32);

VALUE @cs NVARCHAR = ComponentProperty(@table'FieldCoordSystem.' + @field);

--

-- target table

CREATE TABLE [Grid Table] (

    [tileX] INT32[tileY] INT32,

    [pixelX] INT32[pixelY] INT32,

    [imageX] INT32[imageY] INT32,

    [Cell] GEOM,

    [Value] UINT8X3,   -- source data type (restored)

    INDEX [Y_X_x] BTREE ([imageY][imageX]), -- same order as image

    INDEX [Cell_x] RTREE ([Cell]),

    PROPERTY 'FieldCoordSystem.Cell' @cs -- same as image

    );

-- target drawing

CREATE DRAWING [Grid] (

    PROPERTY 'FieldGeom' 'Cell',

    PROPERTY 'Table' '[Grid Table]'

    );

--

INSERT INTO [Grid table]

    (

    [tileX][tileY],

    [pixelX][pixelY]

    [imageX][imageY],

    [Cell],

    [Value]

    )

SELECT

    [tileX][tileY],

    [pixelX][pixelY]

    [imageX][imageY],

    GeomMakeRect(

        VectorMakeX4([imageX][imageY][imageX] + 1, [imageY] + 1)

        ) AS [Cell],

    --[Value] -- FLOAT64XN

    CAST([Value] AS INT32X3-- restore source data type

FROM

    (

    SELECT

        [tileX][tileY],

        [X] AS [pixelX][Y] AS [pixelY]

        [tileX] * @dimX + [X] AS [imageX],

        [tileY] * @dimY + [Y] AS [imageY],

        [Value]

    FROM

        (

        SELECT [X] AS [tileX][Y] AS [tileY],

        SPLIT CALL TileToValuesX3([Tile])

            --> X, Y, Value (always FLOAT64XN)

        FROM @table

        )

    )

-- exclude pixels outside image Rect

--WHERE [imageX] BETWEEN @minX AND @maxX

--AND [imageY] BETWEEN @minY AND @maxY

;

Attachments:
Extract cells and values (int32x3) forum.sql

Manifold User Community Use Agreement Copyright (C) 2007-2017 Manifold Software Limited. All rights reserved.