Subscribe to this thread
Home - General / All posts - M9 Modifying pixels in a raster using a vector features as analysis masks
danb

2,064 post(s)
#15-May-19 05:16

I can’t remember if I have asked this before but have wanted to be able to do it for a while in M9. I have a vector feature which I want to transpose onto an image to modify or extract the overlain pixels.

Thinking of ways to do this, I know there was forum mention of a method similar to M8’s 'no interpolation' method when pasting a vector feature to a grid, but I haven’t seen anything similar in the current range of interpolators. We also now have features for performing zonal type functions collecting data under a vector feature from an overlain raster, but as far as I am aware, not the other way round.

With this in mind, is it currently possible to extract the pixels under a vector feature or interact with them such as setting their value in SQL, or is this a feature request?

Making (irregular) selections in images and modifying the pixel values in some way is a very powerful feature of M8 which I would love to see in M9.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

Dimitri


7,413 post(s)
#15-May-19 05:48

If you don't mean the raster to vector interaction as shown in Example: Transfer DEM Terrain Heights to Areas in a Drawing then that would be a feature request.

tjhb
10,094 post(s)
#15-May-19 08:09

Dan means the other way around: use geometry to look up underlying pixels, then do something with the pixels.

The "something" could be (a) do raster maths, (b) create a new raster from those pixels, (c) modify an existing raster (the source raster or another) only within the scope of those pixels. (Other things?)

I tend to think selecting pixels is a red herring. It is often the most convenient method of doing those things in 8, but is not necessary in 9.

Anyway, all of these things are possible already in 9, just not as easy as they could be.

An exact example might be good, e.g. a Manifold 8 workflow with simple data. Then we can agree on how it can currently be done in 9, and how it might be made easier/more approachable.

danb

2,064 post(s)
#15-May-19 21:16

Thanks Tim. I will put together some simple data for a use case.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

danb

2,064 post(s)
#16-May-19 01:00

Find attached a M8 test project. It contains main three components [VECTOR], [RASTER] and

[1 UPDATE SELECTED PIXELS].

[VECTOR] and [RASTER] are subsets taken from larger extent components in a M9 project.

[1 UPDATE SELECTED PIXELS] is the way I would approach this in M8 i.e. select all objects in [VECTOR], transfer this selection to [RASTER] and then run this query.

The M9 project creates the surface as a subset of a lidar grid including all pixels with heights <=1.8m.and all other pixels set to invisible. The vector feature is the result of a ‘Trace Areas’ transform operating on the same grid with some polygons such as those along the river removed from the resulting drawing.

For this example, all I wish to do is use Manifold 9 to update the pixels under the vector polygons with a value of 99. As Tim suggests however, I am also interested in more general cases of interactions between geom and pixels as outlined.

Attachments:
M8 FORUM.map


Landsystems Ltd ... Know your land | www.landsystems.co.nz

tjhb
10,094 post(s)
#17-May-19 01:13

Cool data!

tjhb
10,094 post(s)
#17-May-19 02:20

Here are the steps I am thinking of, with some effort at optimization.

  • Reproject areas to CS of image
  • Subdivide areas into convex parts, with relevant attribute; result to temporary table T with RTREE on part
  • Make tile boxes for image; result to temporary table U with RTREE on box
  • Filter boxes in U overlapping some part in T
  • Split filtered tiles into pixels, drawing pixel boxes; result to temporary table V with RTREE on box
  • Filter pixel boxes in V overlapping some part in T; result to temporary table W with BTREEDUP on XY tile index
  • Adjust value of joined pixels in W, applying attribute from matched part
  • Left join tiles in U to adjusted pixels in W (join on tile index)
  • Recombine joined tiles with coalesce: adjusted pixel if any, else original
  • Substitute recombined tiles for original tiles as applicable
tgazzard
146 post(s)
#13-Jan-20 00:14

I would like to update the pixel values in a raster using the values in an area vector drawing. Is this now possible using a transform or sql approach in m9?

I have looked through the manual and looked at the latest video on closest rasters which does some interaction with vector points. However, it isn't obvious to me how this would work in m9. This is something that I regularly undertook in m8. Any help appreciated.

tjhb
10,094 post(s)
#13-Jan-20 00:22

It’s not built-in yet Tim.

I think this is Dan’s top outstanding feature request, mine too.

Besides that, it’s definitely time I had a proper go at writing as good a solution as I can.

Overdue. No excuses...

tgazzard
146 post(s)
#13-Jan-20 00:27

Thanks for confirmation that there isn't a feature (yet) that does this.

If you do write a solution, I would love to see this.

dyalsjas
157 post(s)
#13-Jan-20 00:32

Concur, this is a needed capability when working with GIS and images.

Would also like to transform image pixels to points (and back) e.g. elevation image to x,y,z points and the reverse, x,y,z points to elevation image.

Kriging approaches the x,y,z points to elevation image, but direct translation (without interpolation) is sometimes necessary.

Dimitri


7,413 post(s)
#13-Jan-20 04:57

You can do the (and back) part now, transferring pixel values from an image to points in a drawing. See this example.

tjhb
10,094 post(s)
#13-Jan-20 05:30

Yes. That is the reverse, the "and back" part.

The missing part is to transfer data from geometry (of any kind) to coincident pixels.

dyalsjas
157 post(s)
#13-Jan-20 11:14

I reviewed this example, but I'm unable to find an easy way to to create a drawing of points coincident with each pixel in an image. This effort is non-trivial with high resolution DEMs,

Dimitri


7,413 post(s)
#13-Jan-20 12:35

Ah, I see. I misunderstood your interest. I thought that given a raster image and a layer of points above that, you wanted to transfer the value of each pixel under a given point into a field in that point, for example, the way Height values from the raster in the Transform: Closest Raster topic were transferred into the eight points in the Points drawing.

Creating a dense grid of points to overlay a raster is not that difficult. It's tedious, but can be done remarkably quickly given a series of copy/paste/shift and repeat operations, doubling the size of the point swaths that are copied/pasted each time. There have been threads on that.

But, could it be there are easier ways to do what you want than to transform a raster into points?

adamw


10,447 post(s)
#14-Jan-20 16:18

Just wondering, do you need to do this for a couple of points or for massive amounts of points (are we talking about altering some points on a raster or basically re-creating it)?

Both things can be done using a script, but if it is the latter, then maybe there's a different approach you can take and it would be interesting to know why do you want to go from raster to point-for-each-pixel and then back to raster - the first guess is that you want to do something in the middle and cannot find a suitable combination of raster functions to do that in raster (and if so, maybe what we need is a couple more raster functions).

In general, we do have 90% of the machinery to go from vector to raster and will likely do the remaining 10% fairly soon. We have that on the radar for sure.

dyalsjas
157 post(s)
#15-Jan-20 01:49

I am trying to recreate a model that compares surface variability factors across ranges of thirty and three hundred meter diameters. This model can be GPU accelerated when using 10m to 30m elevation models, but I haven’t been able to identify a software that will allow me to develop the same model to work on high resolution elevation models like LiDAR data sets.

CUDA raster acceleration seems to drop off after an analysis mask exceeds a modest radius. This model would have 15 and 150 pixel radius masks respectively when using 1m DEMs.

I thought I might try a brute force approach by converting an elevation model to x,y,x points and using a combination of buffers and queries. The underlying math capability already seems built into Manifold.

I also thought I might explore if a combination of file export and import methods might allow me to get to a suitable set of points.

The final result of the model is effectively a raster image of classes similar to a Land Use Land Cover image. I’m sure I could export a point file from Manifold into something I could kludge back into an image using GDAL or ESRI or something.

I'll once again ask for an option to select units (degrees or percentage) when calculating a slope image.

adamw


10,447 post(s)
#15-Jan-20 08:21

Thanks, that helps.

So, the back and forth between raster and vector is indeed fairly massive, but that's because traditional approaches in the raster seem too slow (150 pixel windows are quite slow, that's true), you want to experiment using geometry.

Let's see what you can do right now:

You can convert a raster into a drawing - that's a simple query. (See TileToValues functions, say if you want an example query.)

You can then analyze that drawing using transforms in 9. Or using a script.

You can also export the drawing to, say, SHP, or perhaps export the table to CSV with X-Y-Z values, and process that data in other software, then import the result back into 9.

After you have the result, if it is just pixels and the pixels are still in their integer positions with just the values changed, it is fairly simple to put them back into a raster using a script. Again, say if you want example code.

If the result is areas or lines, then perhaps wait until we do vector to raster as a built-in tool.

PS: Request for degrees or percentages in slopes noted - I'll bump the wishlist item.

dyalsjas
157 post(s)
#16-Feb-20 00:43

Adam,

It's taken a while for me to get back to this posting.

I've reviewed the various topics here in the user forum that included examples of how to create a table from tile values using the various tiletovalue functions.

I haven't been able to get a working query.

Can you (or any of the other forum members) assist with an example query of how to create a table of x, y, and z values from a single band image. The image I'm currently working with is an image of slope values.

I'm in the middle a trying an x.y.x file export, but it's surprisingly slow.

v/r

jd

adamw


10,447 post(s)
#17-Feb-20 08:52

Sure.

Given an image named 'german_alps', with single-channel tiles of the standard 128x128 size, this will produce the table with x - y - value for each pixel:

--SQL9

VALUE @tilesizex INT32 = 128;

VALUE @tilesizey INT32 = 128;

SELECT tilex * @tilesizex + x AS x, tiley * @tilesizey + y AS y, value FROM (

  SELECT x AS tilex, y AS tiley, SPLIT CALL TileToValues(tile)

  FROM [german_alps]

);

The inner query splits all tiles into pixels and keeps tile X and tile Y for each pixel. The outer query combines tile X / Y and pixel X / Y within a tile into pixel X / Y within an image, and keeps pixel value.

You can convert the outer query into SELECT ... INTO to put this into a table. If you do that, it might be worth it to add MFD_ID / MFD_ID_X to that table to make it editable - or you can create a BTREE index on X / Y (the combinations of which are unique), that will make the table editable as well.

dyalsjas
157 post(s)
#17-Feb-20 19:21

I've tried this query, but I'm not sure the image I'm working with has standard tile sizes.

The tile description shows <tile, 256x16, float32>

When running the query above, the value column in the query result shows NULL.

I've tried updating the tile sizes to 256, but still receive NULL within the value column of the query result.

tjhb
10,094 post(s)
#17-Feb-20 21:06

text

text

When running the query above, the value column in the query result shows NULL.

That might be normal...

(1) How far are you scrolling down? If you are scrolling all the way to the bottom, remember that that means you get the "first" 50,000 records, which here means the "first" 50,000 pixels. ("First" here is arbitrary in theory, but they will usually be drawn from the tiles with the lowest indexes, though this could be disrupted.)

(2) Could the displayed pixels all be <null>, and this still be normal? Yes, if you have a large number of empty tiles. That is, tiles that are present, but contain only invisible pixels. That can be the case with images imported/linked from other formats, including Manifold 8. Manifold 9 is normally more efficient, since it normally only writes tiles which have some visible pixel(s).

To test more definitely, just add a filter to Adam's outer query.

SELECT...

FROM

(

SELECT...

)

WHERE [Value] IS NOT NULL;

That will return up to 50,000 records showing visible pixels... or no records at all, if all pixels were invisible.

adamw


10,447 post(s)
#18-Feb-20 07:58

I've tried updating the tile sizes to 256, but still receive NULL within the value column of the query result.

Maybe that's what you did, but just in case - if your tiles are 256x16, set @tilesizex to 256 and @tilesizey to 16, do not set both to 256.

As Tim says, NULL values represent invisible pixels. If you don't want them, you can filter them out with WHERE.

dyalsjas
157 post(s)
#22-Mar-20 17:55

I've just revisited this query, (including the filter to exclude NULL values), and the output is what I expected.

My next question is: How do I have the X and Y values reflect the projected coordinates of each pixel in the image instead of the pixel coordinates?

adamw


10,447 post(s)
#03-May-20 08:21

A very late reply, I somehow missed the question in March, sorry.

The pixel coordinates adjusted for tile number in the query I posted *are* the coordinates in the coordinate system of the image. They are translated to physical units (eg, meters) using local scale and local offset parameters which you can see at the bottom of the coordinate system dialog for the image. You can translate them to meters / etc, in the query, but you normally don't really need to. Usually what you want to do with them is transform to the coordinate system of a different component - this can be done using coordinate converter objects, and the pixel coordinates are perfectly good for that. (Tell me if you need an example of translating pixel coordinates to the coordinate system of a different component.)

Anyway, I believe a lot of what we have been doing in the thread can now be done easier with the functions and the extension of the Join dialog that we did in:

Manifold System 9.0.171.4

Take a look at the comments, too, they contain various example queries.

dyalsjas
157 post(s)
#08-May-20 00:48

Adam,

I missed your reply earlier this week as well.

Yes please I could use some help with this.

I'm not sure about the schema for a table to SELECT ... INTO

I would love an example of translating the pixel coordinates of the query output to the coordinate system of the source image.

As an example set of data to work with, I'm attaching an extract of the trace_areas_SJC project available for download from Manifold.net.

My very minor changes to your script are included as well.

Attachments:
LULC_Example.mxb

adamw


10,447 post(s)
#08-May-20 09:01

Here you go:

--SQL9

 

VALUE @conv TABLE = CALL CoordConverterMake(

  ComponentCoordSystemAutoXY([Drawing]),       -- target

  ComponentCoordSystemAutoXY([grid_cell LULC]-- source

);

 

VALUE @tilesizex INT32 = 128;

VALUE @tilesizey INT32 = 128;

 

SELECT CoordConvertPoint(@conv,

  VectorMakeX2(tilex * @tilesizex + x, tiley * @tilesizey + y)) AS xy,

  value

FROM (

  SELECT x AS tilex, y AS tiley, SPLIT CALL TileToValues(Tile)

  FROM [grid_cell LULC]

WHERE value IS NOT NULL;

We take the X and Y pixel coordinates, create a point, then convert that point from the coordinate system of the image to the coordinate system of a different component (I added a drawing as an illustration).

Extended MXB attached.

Attachments:
LULC_Example-2.mxb

dyalsjas
157 post(s)
#08-May-20 17:55

Adam,

I'm attaching a further update to our exchanges. To keep it distinct, I'll start naming it by date (and exchange number if necessary).

When I run your query, above, in the .mxb you shared, (with a small mod to output the query to a table), the resulting table shows a coordinates attribute with lon / lat as a decimal degree pair and a value attribute.

If I reproject the target drawing to the coordinate system of the image (UTM 10N), the output table shows the pixel coordinate, not the UTM 10N coordinates.

I suspect I'm misunderstanding something fundamental in how Manifold deals with coordinate systems.

Two further questions:

Is the VectorMakeX2 function required as part of the CoordConverterMake process or can the converted X and Y coordinates be output to separate attributes as in the first query you shared?

Rather than calling a projection from another drawing, could a projection be defined as a VALUE that can be called later in the query?

Attachments:
2020_0508_LULC_Example_1.mxb

adamw


10,447 post(s)
#09-May-20 08:20

The coordinate system of the image is UTM + scales and offsets designed so that the lower left corner of the lower left pixel is at 0, 0, the lower left corner of the next pixel to the right is at 1, 0, etc. You can see the scales and offsets applied to UTM in the coordinate system dialog, in the metrics section at the bottom.

If you copy the coordinate system of the image to the drawing, the coordinate system of the drawing becomes the same, with 0, 0 being where the lower left corner of the lower left pixel of the image is, etc. It's not just UTM, zone X, it's UTM, zone X + adjustments. If you want to see coordinates in UTM, zone X, without adjustments, you should project to the coordinate system which does not have those adjustments. One way to do this is to edit the coordinate system of the drawing, edit metrics and set local scales to 1, 1 and local offsets to 0, 0 (all in meters). Alternatively, you can do this dynamically in a query. I altered the 'TileToValues Project to UTM Drawing' in the attached file to do the latter.

The x2 value returned by the CoordConvertPoint call can, of course, be split into X and Y components. To do this, use VectorValue. I altered both the 'TileToValues Project to UTM Drawing' and the 'TileToValues Project to Geo Drawing' queries to do this during SELECT INTO.

A projection can be defined as a value, absolutely. Eg:

--SQL9

VALUE @target NVARCHAR = ComponentCoordSystemAutoXY([component]);

  -- or

VALUE @target NVARCHAR = 'EPSG:3857';

New file attached.

Attachments:
2020_0509_LULC_Example_1.mxb

danb

2,064 post(s)
#15-Jan-20 04:36

I think this is Dan’s top outstanding feature request, mine too.

Right on the money Tim. Definitely my number one currently.

I have been away for a while but what a treat to come back to. Two great builds.

I am very much looking forward to really getting to grips with the distance tools and the las library looks like a game changer for me and the way we tend to store las data (primarily to cater to legacy GIS).

The improvements to day to day usability such as 'rename related' are just great.


Landsystems Ltd ... Know your land | www.landsystems.co.nz

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