[gdal-dev] Questions about polygonization
I am trying to use GDAL's polygonize algorithms to help me identify regions in landcover data. I made a trivial example, but I am having trouble understanding the results that I get, nor can I determine how to get the results I want. Given the following raster: 11 11 11 11 11 11 11 11 11 11 11 11 12 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 If I were to polygonize this raster, I would expect two polygons: one which contains the single point 12, and another with an outer ring was the perimeter of the raster and an inner ring circling the single point 12. I would also expect that the first polygon would contain the single point 12 and that the second polygon would contain all the points 11. Instead, what I got was a polygon starting from the single point 12 and extending down and to the right one pixel, and another polygon with an outer ring that went off the down and right edges and an inner ring equal to the first polygon's outer ring. The first polygon doesn't actually contain any points and the second one doesn't contain any points on either ring but does contain all the other points. What do I need to do to get the results I am expecting? Jack. -- mathuin at gmail dot com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Questions about polygonization
Hi Jack,
I'm not sure if I understand the description of the shapes you see.
Can you provide the WKT for them? I seem to get two polygons that look
like they should:
import numpy as np
import rasterio.features
from shapely.geometry import shape
ar = np.ones((6, 5), 'B') * 11
ar[2, 2] = 12
gt = [0, 1, 0, 6, 0, -1]
for (sh, v) in rasterio.features.shapes(ar, transform=gt):
print('%s: %s' % (v, shape(sh)))
11: POLYGON ((0 6, 0 0, 5 0, 5 6, 0 6), (2 4, 3 4, 3 3, 2 3, 2 4))
12: POLYGON ((2 4, 2 3, 3 3, 3 4, 2 4))
-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Questions about polygonization
Jack, The GDALPolygonize algorithm draws the polygon edges along the pixel edges. It assumes the pixels to be rectangular areas instead of points. In the example you described, the first polygon should contain four vertices with the pixel in the centre. Make sure that the program you are using to display them is not shifting the pixels. On 20 Aug 2014 07:12, "John Twilley" wrote: > I am trying to use GDAL's polygonize algorithms to help me identify > regions in landcover data. I made a trivial example, but I am having > trouble understanding the results that I get, nor can I determine how > to get the results I want. > > Given the following raster: > > 11 11 11 11 11 > 11 11 11 11 11 > 11 11 12 11 11 > 11 11 11 11 11 > 11 11 11 11 11 > 11 11 11 11 11 > > If I were to polygonize this raster, I would expect two polygons: one > which contains the single point 12, and another with an outer ring was > the perimeter of the raster and an inner ring circling the single > point 12. I would also expect that the first polygon would contain > the single point 12 and that the second polygon would contain all the > points 11. > > Instead, what I got was a polygon starting from the single point 12 > and extending down and to the right one pixel, and another polygon > with an outer ring that went off the down and right edges and an inner > ring equal to the first polygon's outer ring. The first polygon > doesn't actually contain any points and the second one doesn't contain > any points on either ring but does contain all the other points. > > What do I need to do to get the results I am expecting? > > Jack. > -- > mathuin at gmail dot com > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Questions about polygonization
Hi Mike!
I get the same polygons you get -- what I don't understand is why I'm
getting those polygons instead of the ones that I expected. Can you
explain why those polygons look like they should?
Jack.
--
mathuin at gmail dot com
On Tue, Aug 19, 2014 at 7:15 PM, Mike Toews wrote:
> Hi Jack,
>
> I'm not sure if I understand the description of the shapes you see.
> Can you provide the WKT for them? I seem to get two polygons that look
> like they should:
>
> import numpy as np
> import rasterio.features
> from shapely.geometry import shape
> ar = np.ones((6, 5), 'B') * 11
> ar[2, 2] = 12
> gt = [0, 1, 0, 6, 0, -1]
> for (sh, v) in rasterio.features.shapes(ar, transform=gt):
> print('%s: %s' % (v, shape(sh)))
>
> 11: POLYGON ((0 6, 0 0, 5 0, 5 6, 0 6), (2 4, 3 4, 3 3, 2 3, 2 4))
> 12: POLYGON ((2 4, 2 3, 3 3, 3 4, 2 4))
>
> -Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Questions about polygonization
The program is not shifting the pixels, but the program does assume the pixels to be points and not rectangular areas so maybe I will have to shift the pixels. Do you know what transformation I would need to use? Jack. -- mathuin at gmail dot com On Tue, Aug 19, 2014 at 7:17 PM, Chaitanya kumar CH wrote: > Jack, > The GDALPolygonize algorithm draws the polygon edges along the pixel edges. > It assumes the pixels to be rectangular areas instead of points. > In the example you described, the first polygon should contain four vertices > with the pixel in the centre. > Make sure that the program you are using to display them is not shifting the > pixels. > > On 20 Aug 2014 07:12, "John Twilley" wrote: >> >> I am trying to use GDAL's polygonize algorithms to help me identify >> regions in landcover data. I made a trivial example, but I am having >> trouble understanding the results that I get, nor can I determine how >> to get the results I want. >> >> Given the following raster: >> >> 11 11 11 11 11 >> 11 11 11 11 11 >> 11 11 12 11 11 >> 11 11 11 11 11 >> 11 11 11 11 11 >> 11 11 11 11 11 >> >> If I were to polygonize this raster, I would expect two polygons: one >> which contains the single point 12, and another with an outer ring was >> the perimeter of the raster and an inner ring circling the single >> point 12. I would also expect that the first polygon would contain >> the single point 12 and that the second polygon would contain all the >> points 11. >> >> Instead, what I got was a polygon starting from the single point 12 >> and extending down and to the right one pixel, and another polygon >> with an outer ring that went off the down and right edges and an inner >> ring equal to the first polygon's outer ring. The first polygon >> doesn't actually contain any points and the second one doesn't contain >> any points on either ring but does contain all the other points. >> >> What do I need to do to get the results I am expecting? >> >> Jack. >> -- >> mathuin at gmail dot com >> ___ >> gdal-dev mailing list >> gdal-dev@lists.osgeo.org >> http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Questions about polygonization
Try qgis to display the raster and the vector output. It shows the raster pixels as rectangular areas. To confirm that there is no shift, compare their extents using gdalinfo and ogrinfo. On 20 Aug 2014 07:56, "John Twilley" wrote: > The program is not shifting the pixels, but the program does assume > the pixels to be points and not rectangular areas so maybe I will have > to shift the pixels. Do you know what transformation I would need to > use? > > Jack. > -- > mathuin at gmail dot com > > > On Tue, Aug 19, 2014 at 7:17 PM, Chaitanya kumar CH > wrote: > > Jack, > > The GDALPolygonize algorithm draws the polygon edges along the pixel > edges. > > It assumes the pixels to be rectangular areas instead of points. > > In the example you described, the first polygon should contain four > vertices > > with the pixel in the centre. > > Make sure that the program you are using to display them is not shifting > the > > pixels. > > > > On 20 Aug 2014 07:12, "John Twilley" wrote: > >> > >> I am trying to use GDAL's polygonize algorithms to help me identify > >> regions in landcover data. I made a trivial example, but I am having > >> trouble understanding the results that I get, nor can I determine how > >> to get the results I want. > >> > >> Given the following raster: > >> > >> 11 11 11 11 11 > >> 11 11 11 11 11 > >> 11 11 12 11 11 > >> 11 11 11 11 11 > >> 11 11 11 11 11 > >> 11 11 11 11 11 > >> > >> If I were to polygonize this raster, I would expect two polygons: one > >> which contains the single point 12, and another with an outer ring was > >> the perimeter of the raster and an inner ring circling the single > >> point 12. I would also expect that the first polygon would contain > >> the single point 12 and that the second polygon would contain all the > >> points 11. > >> > >> Instead, what I got was a polygon starting from the single point 12 > >> and extending down and to the right one pixel, and another polygon > >> with an outer ring that went off the down and right edges and an inner > >> ring equal to the first polygon's outer ring. The first polygon > >> doesn't actually contain any points and the second one doesn't contain > >> any points on either ring but does contain all the other points. > >> > >> What do I need to do to get the results I am expecting? > >> > >> Jack. > >> -- > >> mathuin at gmail dot com > >> ___ > >> gdal-dev mailing list > >> gdal-dev@lists.osgeo.org > >> http://lists.osgeo.org/mailman/listinfo/gdal-dev > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Questions about polygonization
Using qgis helped me see what you are saying. The raster pixels are being treated as rectangular areas but the actual coordinates are the *corners* of those areas. Modifying my code to offset the coordinates by 0.5 causes the Contains test (and others) to produce results that meet my expectations. The really frustrating thing is that I think I had to learn this already at least once over the past decade or so I've been tinkering with GDAL. I can't find a particularly embarrassing post in the archives, but I'm confident it's there. :-) Thank you for your patience. -- mathuin at gmail dot com On Tue, Aug 19, 2014 at 7:32 PM, Chaitanya kumar CH wrote: > Try qgis to display the raster and the vector output. It shows the raster > pixels as rectangular areas. To confirm that there is no shift, compare > their extents using gdalinfo and ogrinfo. > > On 20 Aug 2014 07:56, "John Twilley" wrote: >> >> The program is not shifting the pixels, but the program does assume >> the pixels to be points and not rectangular areas so maybe I will have >> to shift the pixels. Do you know what transformation I would need to >> use? >> >> Jack. >> -- >> mathuin at gmail dot com >> >> >> On Tue, Aug 19, 2014 at 7:17 PM, Chaitanya kumar CH >> wrote: >> > Jack, >> > The GDALPolygonize algorithm draws the polygon edges along the pixel >> > edges. >> > It assumes the pixels to be rectangular areas instead of points. >> > In the example you described, the first polygon should contain four >> > vertices >> > with the pixel in the centre. >> > Make sure that the program you are using to display them is not shifting >> > the >> > pixels. >> > >> > On 20 Aug 2014 07:12, "John Twilley" wrote: >> >> >> >> I am trying to use GDAL's polygonize algorithms to help me identify >> >> regions in landcover data. I made a trivial example, but I am having >> >> trouble understanding the results that I get, nor can I determine how >> >> to get the results I want. >> >> >> >> Given the following raster: >> >> >> >> 11 11 11 11 11 >> >> 11 11 11 11 11 >> >> 11 11 12 11 11 >> >> 11 11 11 11 11 >> >> 11 11 11 11 11 >> >> 11 11 11 11 11 >> >> >> >> If I were to polygonize this raster, I would expect two polygons: one >> >> which contains the single point 12, and another with an outer ring was >> >> the perimeter of the raster and an inner ring circling the single >> >> point 12. I would also expect that the first polygon would contain >> >> the single point 12 and that the second polygon would contain all the >> >> points 11. >> >> >> >> Instead, what I got was a polygon starting from the single point 12 >> >> and extending down and to the right one pixel, and another polygon >> >> with an outer ring that went off the down and right edges and an inner >> >> ring equal to the first polygon's outer ring. The first polygon >> >> doesn't actually contain any points and the second one doesn't contain >> >> any points on either ring but does contain all the other points. >> >> >> >> What do I need to do to get the results I am expecting? >> >> >> >> Jack. >> >> -- >> >> mathuin at gmail dot com >> >> ___ >> >> gdal-dev mailing list >> >> gdal-dev@lists.osgeo.org >> >> http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
