0

I am hoping for advice as to how I can convert GeoJSON geometry into a tif image using rasterio. I have tried a lot of things, but all of them do not rasterize all the shapes found in the GeoJSON (more like 80% of the file is rasterized). How can I ensure all the geometry is rasterized and is of adequate size? Let me know if my question is unclear.

DarkDrassher34
  • 69
  • 3
  • 15
  • Hi, Could you please share, if possible, an example of a GeoJSON that you would like to rasterize, and the code you have written so far that only rasterizes 80% of the shapes? This could give us a starting point. – guampi Dec 10 '19 at 13:31
  • I cannot share the GeoJSON, but it is just one big `MultiPolygon` object. I have `gdal` installed, so the code I used was: `gdal_rasterize -b 255 -b 0 -b 0 -at -ts 4450 4900`. Hope this helps. – DarkDrassher34 Dec 10 '19 at 13:42
  • I guess there's a typo in your previous comment, and you meant something like `-burn 255` instead of `-b 255`. Apart from the typo, I don't see any reason why the command you posted would not work. Is your question how to make it work with GDAL, or how you can translate this into rasterio code? – guampi Dec 10 '19 at 13:51
  • Well, it works. It rasterizes the `MultiPolygon` feature, but the problem I think is with the parameter (`ts`). It is arbitrarily picked. I need to know how I can choose the right parameters to then run `gdal_rasterize`. As it stands, not all the geometry is being burned. – DarkDrassher34 Dec 10 '19 at 13:52
  • I'd be curious to see if you can provide a toy example that showcases your error because I don't manage to reproduce a case where only part of the geometry is being burned. – guampi Dec 10 '19 at 14:05
  • Actually, maybe your problem stems from the fact that you choose a raster size, through `-ts`, that is not resolute enough to show all of your features. See for example this MultiPolygon: https://gist.github.com/glostis/e6a6bab1f114cc99222a98d3bc10f419 If you run `gdal_rasterize -burn 255 -ts 100 100 map.geojson a.tif` you have enough resolution with 100 pixels to see the two features, but if you run `gdal_rasterize -burn 255 -ts 10 10 map.geojson a.tif` the image resolution is too low, and the lower left feature looks strange due to that. – guampi Dec 10 '19 at 14:14
  • Exactly. That is what I am getting at. How do I choose that parameter? I would think `gdal` would do some simple algorithmic calculation, but it does not. – DarkDrassher34 Dec 10 '19 at 14:19
  • You are rasterizing vectors here, so IMO `gdal` cannot choose the resolution for you, it is up to you to determine what raster resolution you want. FYI, you might be interested in the `-tr` option with which you directly set the resolution you want, instead of `-ts` with which you set the size. – guampi Dec 10 '19 at 14:23
  • I have also tried it. I think my problem here is that my `MultiPolygon` object is also in `lat/lon` coordinates, which makes it harder for the resolution to be 1:1. I know this sounds idealistic, but how do you go about setting that resolution? If there is not an answer other than 'keep trying' then I think I have stumbled into a research problem. – DarkDrassher34 Dec 10 '19 at 14:29
  • If you do not want to lose too many details in the rasterization process, I guess a good rule of thumb would be to set the resolution to be lower than the typical size of your individual features. For example, if your features look like squares ~1km large, a "good" resolution would be in the 10-100m range. A way to programmatically get the typical size of your features would be to compute their `minimum_rotated_rectangle` for example: https://shapely.readthedocs.io/en/stable/manual.html#object.minimum_rotated_rectangle – guampi Dec 10 '19 at 14:46
  • Great. Copy your last comment in the answer, so I can give you some points. – DarkDrassher34 Dec 10 '19 at 14:53

1 Answers1

1

Your problem stems from the fact that when you rasterize a shape, using gdal_translate for example, you must determine the resolution of the raster, and it must be chosen in accordance to the "size" of your feature vectors if you want your raster to retain enough information.

If you do not want to lose too many details in the rasterization process, I guess a good rule of thumb would be to set the resolution to be lower than the typical size of your individual features. For example, if your features look like squares ~1km large, a "good" resolution would be in the 10-100m range. A way to programmatically get the typical size of your features would be to compute their minimum_rotated_rectangle for example: https://shapely.readthedocs.io/en/stable/manual.html#object.minimum_rotated_rectangle

guampi
  • 306
  • 1
  • 8
  • Ok. You calculate the `minimum_rotated_rectangle` and then get the area from that? Or what gives me the parameter? – DarkDrassher34 Dec 10 '19 at 16:04
  • You can take the size of the smallest side of your rectangle, for example, to get a feel for how small your shapes can be. Then take a resolution smaller than that length, in order for your smallest shape to be larger than just one pixel in your raster. Once again this is just a rule of thumb, it should not be taken as the golden truth. – guampi Dec 10 '19 at 16:12
  • A rule of thumb is usually better than being in the dark. I will investigate further. Thanks. – DarkDrassher34 Dec 10 '19 at 16:17