I have a TIFF file and I'm trying to use rasterstats to get statistics of different geometries defined in a shapefile or a TopoJSON file. But I'm getting the count as 0 for all the geometries. My guess is that it's because the shapefile doesn't align properly with the raster TIFF file, maybe I'm wrong.

Below is the output of rio info:

    "tiled": true,
    "driver": "GTiff",
    "width": 3390,
    "nodata": 32767.0,
    "colorinterp": ["grey"],
    "lnglat": [94.65791519501956, 28.06018778470516],
    "interleave": "band",
    "crs": "+a=6371007.181 +b=6371007.181 +lon_0=0 +no_defs +proj=sinu +units=m +x_0=0 +y_0=0",
    "blockxsize": 128,
    "dtype": "int16",
    "shape": [1348, 3390],
    "height": 1348,
    "count": 1,
    "blockysize": 128,
    "res": [231.65635826395808, 231.65635826395797],
    "bounds": [8895604.157333, 2964017.6537250075, 9680919.211847818, 3276290.424664823],
    "transform": [231.65635826395808, 0.0, 8895604.157333, 0.0, -231.65635826395797, 3276290.424664823]

And below are some of the relevant fields from the TopoJSON file:

    "type": "Topology",
    "objects": {
        "ArunachalNewDistrict-wise": {
            "type": "GeometryCollection",
            "bbox": [91.55482482910162, 26.656021118164176, 97.41516113281256, 29.463344573974723],
        "geometries": [{...}, {...}, ...],
        "arcs": [...],
        "transform": {
        "scale": [0.0005860922395950532, 0.000280760421623217],
        "translate": [91.55482482910162, 26.656021118164176]

UPDATE: Here's the fio info of the shapefile also:

    "schema": {
        "geometry": "Polygon",
        "properties": {
            "ID_0": "int:10",
            "ISO": "str:3",
            "NAME_0": "str:75",
            "ID_1": "int:10",
            "NAME_1": "str:75",
            "ID_2": "int:10",
            "NAME_2": "str:75",
            "ID_3": "int:10",
            "NAME_3": "str:75",
            "TYPE_3": "str:50",
            "ENGTYPE_3": "str:50",
            "NL_NAME_3": "str:75",
            "VARNAME_3": "str:100"
    "driver": "ESRI Shapefile",
    "crs_wkt": "GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295],AUTHORITY[\"EPSG\",\"4326\"]]",
    "name": "ArunachalNewDistrict-wise",
    "crs": "EPSG:4326",
    "bounds": [91.55482482910162, 26.656021118164176, 97.41516113281256, 29.463344573974723],
    "count": 15

I'm completely new to this, so most probably I'm missing something obvious. Can anybody tell me what's going on?

I can provide more info if necessary.

  • The only way for knowing if both layers are aligned properly is putting them together at the Map Canvas. – xunilk Jun 17 '17 at 16:24

The only way for knowing if both layers are aligned properly is putting them together at the Map Canvas. I did that (see next image), installed rasterstats python module with easy_install and, finally, ran a similar example as in your link.

>>>from rasterstats import zonal_stats
>>>stats = zonal_stats("/home/zeito/pyqgis_data/new_polygon.shp", "/home/zeito/pyqgis_data/utah_demUTM2.tif")
[{'count': 13004, 'max': 2738.0, 'mean': 2022.9993848046754, 'min': 1520.0}] 

It worked in my case; as it can be watching in next image.

enter image description here

  • Thanks a lot, I tried your suggestion, and in fact, both the layers are very far apart (I had to zoom out a lot). Is there a way to align both of them to each other? Also, the vector layer is very small compared to the raster one, so can I scale it to appropriate proportion? – Devesh Khandelwal Jun 17 '17 at 17:06
  • Nevermind, found it. Will try and update here. Thanks. – Devesh Khandelwal Jun 17 '17 at 17:16

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.