pythonraster
I'm a beginner with python. I would like to extract raster (.tif) values from a set of specific points that I have in a shapefile.
I would like to have the same result as the QGIS Point Sampling Tool plugin.
I can't install on Python libraries like qgis.core arcpy pygrass, but I have available GeoPandas GDAL Fiona Shapely rasterio and others.
How can I do this?
I added that recipe to the rasterio documentation. Since it was such a simple shape, in this case I just unzipped the coords in the single record contained by the shapefile. That is, x, y = zip(*features[0]['coordinates'][0]), and then just plot.
More generally, I use PolygonPatch from descartes, and matplotlib.collections.
import fiona
import rasterio
import rasterio.plot
import matplotlib as mpl
from descartes import PolygonPatch
src = rasterio.open("tests/data/RGB.byte.tif")
with fiona.open("tests/data/box.shp", "r") as shapefile:
features = [feature["geometry"] for feature in shapefile]
rasterio.plot.show((src, 1))
ax = mpl.pyplot.gca()
patches = [PolygonPatch(feature) for feature in features]
ax.add_collection(mpl.collections.PatchCollection(patches))
The appearance of the shapes can be customized with keywords like edgecolor or facecolor passed to PolygonPatch. To produce a thick red line as in the example, replace the last two lines in the example above with:
patches = [PolygonPatch(feature, edgecolor="red", facecolor="none", linewidth=2) for feature in features]
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))
The match_original keyword is necessary in the second example because parameters like facecolor and edgecolor can also be set in PatchCollection, and the default is to ignore settings of the passed in patches in favor of those set by PatchCollection. Doing so using PolygonPatch gives more flexibility to set different colors, widths, etc., for each patch that you add.
You're almost there. Updated docs are at https://fiona.readthedocs.io
Your schema is just a simple dict with 'geometry' and 'properties' (for the attribute fields) items.
You then need to populate a record with the coordinates as x,y tuples (no shapely required for something simple like this) and whatever attributes you like.
The coordinate tuples/lists for a polygon need to be nested:
So a geometry would look something like:
[((x,y),(x,y),(x,y),(x,y))] # Note use of list for single element, can use a tuple, but that requires a trailing comma (((x,y),(x,y etc)),)
Just be careful with coordinate ordering as I've noted in the snippet below.
from collections import OrderedDict
import fiona
from fiona.crs import from_epsg
# example of a single rectangle that we get from that analysis
geolocated_xs = (-113.21645410141903, -113.21645410141903, -113.2167506315605, -113.2167506315605)
geolocated_ys = (26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)
# You could just zip the geolocated_xys together to form your polygon,
# i.e. polygon = [zip(geolocated_xs, geolocated_ys)]
# but you need to be careful with ordering as the ordering of your
# coordinates will give a bow-tie self intersection
# (they are currently in lr->ur->ll->ul order)
# Fix the ordering otherwise you'll get a bow-tie self intersection
lr, ur, ll, ul = zip(geolocated_xs, geolocated_ys)
polygon = [(lr, ur, ul, ll)]
crs = from_epsg(4326)
# Define your schema as a polygon geom with a couple of fields
schema = {
'geometry': 'Polygon',
'properties': OrderedDict([
('somecol', 'str'),
('anothercol', 'int')
])
}
with fiona.open(
'oo.shp',
'w',
driver='ESRI Shapefile',
crs=crs,
schema=schema) as c:
record = {
'geometry': {'coordinates': polygon, 'type': 'Polygon'},
'id': '0',
'properties': OrderedDict([('somecol', 'Some Value'),
('anothercol', 12345)
]),
'type': 'Feature'}
c.write(record)
And for multiple records:
# example of a multiple rectangles using your coordinate ordering
geolocated_xys = (
(
(-113.21645410141903, -113.21645410141903, -113.2167506315605, -113.2167506315605),
(26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)
),
(
(-113.21615757127756, -113.21615757127756, -113.21645410141903, -113.21645410141903),
(26.827766731681987, 26.828098563030775, 26.827766731681987, 26.828098563030775)
)
)
etc...
with fiona.open(
'c:/temp/oo.shp',
'w',
driver='ESRI Shapefile',
crs=crs,
schema=schema) as c:
for id, polygon in enumerate(geolocated_xys):
# Fix the ordering otherwise you'll get a bow-tie self intersection
lr, ur, ll, ul = zip(*polygon)
polygon = [(lr, ur, ul, ll)]
record = {
'geometry': {'coordinates': polygon, 'type': 'Polygon'},
'id': id,
'properties': OrderedDict([('somecol', 'Some Value'),
('anothercol', 12345)
]),
'type': 'Feature'}
c.write(record)
Best Answer
I use rasterio and geopandas. My example uses UTM coordinates. Obviously these fields will depend on your particular shapefile. In my experience this produces indentical results to the QGIS Point Sampling Tool. I like this method because the resulting DataFrame of point and corresponding raster values is easy to analyze (e.g. compute the difference between the point and raster values) and then plot or export to some other tabular format (e.g. CSV).