1

So, I have a TIF georeferenced map and a rather messy shape that I need to digitize. It is green but has a lot of holes in it that are light green color. So, I want to capture all green pixels and draw a polygon. So far, I managed to get all of the green pixels rows and cols as tuples as well as coordinates. I also manage to get polygons based on them but they don't match what is needed. The polygons are basically stacked on the edge of the raster and extending way below it. Click to see image

I read the raster using rasterio and got color pixels using:

rasterio_image = rasterio.open('repro.TIF')

image = cv2.imread('Name.TIF') image =
cv2.cvtColor(image, cv2.COLOR_BGR2RGB)

def RGB2HEX(color):
    return "#{:02x}{:02x}{:02x}".format(int(color[0]), int(color[1]), int(color[2]))

def get_image(image_path):
    image = cv2.imread(image_path)
    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
    return image

def get_colors(image, number_of_colors, show_chart):

    modified_image = cv2.resize(image, (600, 400), interpolation = cv2.INTER_AREA)
    modified_image = modified_image.reshape(modified_image.shape[0]*modified_image.shape[1], 3)

    clf = KMeans(n_clusters = number_of_colors)
    labels = clf.fit_predict(modified_image)

    counts = Counter(labels)
    # sort to ensure correct color percentage
    counts = dict(sorted(counts.items()))

    center_colors = clf.cluster_centers_
    # We get ordered colors by iterating through the keys
    ordered_colors = [center_colors[i] for i in counts.keys()]
    hex_colors = [RGB2HEX(ordered_colors[i]) for i in counts.keys()]
    rgb_colors = [ordered_colors[i] for i in counts.keys()]

    if (show_chart):
        plt.figure(figsize = (8, 6))
        plt.pie(counts.values(), labels = hex_colors, colors = hex_colors)

    return rgb_colors

get_colors(get_image('Name.TIF'), 8, True)

h = input('Input hex: ').lstrip('#') 
RGB_Color = tuple(int(h[i:i+2], 16) for i in (0, 2, 4)) 

indices = np.where(image == RGB_Color) coordinates = zip(indices[0],indices[1]) 
unique_coordinates = list(set(list(coordinates))) 
x = tuple([x[0] for x in unique_coordinates]) 
y = tuple([y[1] for y in unique_coordinates])

some = rasterio.transform.xy(rasterio_image.transform, rows=x, cols=y, 
    offset='center')
wgs_coords = list(zip(some[0], some[1]))
# unique_pixels are coordinates, unique_pixels_true are pixels
unique_pixels = np.array(wgs_coords)
unique_pixels_true = np.array(unique_coordinates)

vectors = 
rasterio.features.shapes(source=unique_pixels.astype(np.float32), 
    transform = rasterio_image.transform)
vectors = list(vectors)
values = [value for polygon, value in vectors]
polygons = [shape(polygon) for polygon, value in vectors]

# Create a geopandas dataframe populated with the polygon shapes
s2_poly_gdf = gpd.GeoDataFrame(data={"id": values}, geometry=polygons, 
    crs=rasterio_image.crs)

s2_poly_gdf.to_file("Test_Polygons.shp")

Why do the polygons appear stacked starting with the left upper corner of the raster?

I think the problem is in x and y since indices has three tuples, one for x, one for y and one for presumably RGB as when printing indices provides:

(array([ 32,  33,  33, ..., 417, 418, 418], dtype=int64), array([242, 
144, 145, ...,  22,  30, 608], dtype=int64), array([1, 1, 1, ..., 0, 1, 
1], dtype=int64))

So, I am probably missing a dimension which is causing the problem later on. Reading with rasterio gets me (431, 636) and with cv2.imread (431, 636, 3). I know there are some guides to vectorizing the entire raster but I want to only vectorize pixels of specific color I identified.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.