2

I am beginner to work with GIS, so be patient.

My first task is to do:

  • Read a raster file (Done)
  • Create a shape file (work on it)
  • Create a new raster from the shapefile (not started)

I am having trouble with to plot the shapefile that I have created over the raster image. The bounding box do not appear at image. How can I figure this out?

I was trying to do something similar to this post: Plot shapefile on top of raster using plot and imshow from matplotlib

from collections import OrderedDict
import rasterio as rio # funcao similar ao R do raster
import rasterio.plot as show
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.collections as cplt
import fiona
import shapely
from descartes import PolygonPatch
from osgeo import gdal

Read the raster file


filename = "C:/Users/adolf/Desktop/Annia_teste/ALF_16JUN26141551_M2AS_R02C1_pansharpen.tif"
dataset = rio.open(filename)


Get the band, width and height from raster data


nband   = dataset.count
nwidth  = dataset.width
nheight = dataset.height


O tipo de dado fica na propriedade .dtypes


=============================================================================


Dataset georeferencing


Item 1 -  Obtendo o valor corresponde aos vertices da imagem em metros


bounds = dataset.bounds    # Get the spatial geopoint
matrix = dataset.transform # Get the affine matriz to relate the pixel with spatial position
crs    = dataset.crs       # Get Coordinate reference system


=============================================================================


Using Fiona to create the shapefile


Define your schema as a polygon geom with a couple of fields


schema = {
    'geometry': 'Polygon',
    'properties': OrderedDict([
        ('somecol', 'str'),
        ('anothercol', 'int')
  ])
}


Create a set of shapefile which are construct by defining a 256 x 256 pixel


square_size = 256
nrow = int(nheight/square_size)
ncol = int(nwidth/square_size)


nrow = 1 #bypass test
ncol = 1


with fiona.open('oo.shp', 'w', driver='ESRI Shapefile',  crs=crs, schema=schema) as c:


id = 0
for col_id in range(0,ncol):
    for row_id in range(0,nrow):

        #Criando o poligo a partir dos vertices
        right_bottom = matrix*(row_id*square_size,col_id*square_size)
        left_bottom  = matrix*(row_id*square_size,(col_id+1)*square_size)
        right_top    = matrix*((row_id+1)*square_size,col_id*square_size)
        left_top     = matrix*((row_id+1)*square_size,(col_id+1)*square_size)

        # Corrigindo a ordem p/ escrever os pontos na sequencia correta
        #lr, ur, ll, ul = zip(*polygon)
        polygon = [(right_bottom, right_top, left_top, left_bottom)]

        record = {
            'geometry': {'coordinates': polygon, 'type': 'Polygon'},
            'id': id,
            'properties': OrderedDict([('somecol', 'Some Value'),
                                       ('anothercol', 12345)
                                       ]),
            'type': 'Feature'}
        #Escrevendo no artigo
        c.write(record)
        id = id + 1




Reading the shape file I have just created


with fiona.open("oo.shp", "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]


patches = [PolygonPatch(feature) for feature in features]


ax = rio.plot.show(dataset, adjust='linear')
ax.add_collection(cplt.PatchCollection(patches))
Vince
  • 20,017
  • 15
  • 45
  • 64
Adolfo Graciano Da Silva
  • 21
  • 2
  • There's a lot of code there to pick an error in. I'm not a shapely or fiona user but I can see you're not closing the polygon, you need to specify your origin point again as the last point to close it: polygon = [(right_bottom, right_top, left_top, left_bottom,right_bottom)], though the order of vertices is counter clockwise which may also be an issue. I'm not sure about the matrix part, I think it would work just as well without matrix, as the example https://fiona.readthedocs.io/en/latest/manual.html#record-geometry the coordinates are simple tuples. – Michael Stimson Aug 17 '20 at 04:05

0 Answers0