I'm new to GIS.

I'm a physicist and I have a simulation which I usually run in numpy arrays (usually ending up with a heatmap or contour plot like below).

enter image description here

However, I now want to run the simulation over the geography of a country or region rather than just a square. I have the shape file for the region I'm interested in, and I know how to plot basic lat/long points, but how can I create some sort of array or grid to run a simulation over?

I've tried turning the shapefile into a binary mask and putting that into a numpy array, but that gets rid of any spatial information.

I'm solving a differential equation over a discrete scalar field (i.e. my numpy array is a discrete grid over some region, and the values in the array are evolved over time). I'm interested in how the geometry of the system affects the answer. Originally I used a square array and experimented with simple binary masks like circles etc., but the model is broadly to do with population movement so I want to use the geometry of a real region to investigate things further.

  • You only need to rasterize shapefile with a dimension equivalent of your array (it will define cell size). You also get 'a square' but outside shapefile values are nodata and inside it the values obtained from simulation. Afterward, you need to convert array in a raster for incorporating spatial information with gdal python module. – xunilk Jan 24 at 18:02
  • Can you tell something more about the process? Seems like an interesting issue to solve but I'm not sure what exactly you want to achieve. Maybe raster algebra is what you're looking for? In that case rasterizing is the first step, just as @xunilk suggested. – Leon Powałka Jan 24 at 21:24
  • I would also go for a raster-based solution, as simulation usually use finite elements. – s.k Jan 24 at 21:27
  • Hi, I've edited with more details about what I'm doing. Would you be able to point me in the right direction for rastering in this context? I've seen generic "how to raster a shapefile" posts here but I don't really understand how to get the kind of "masked-array-but-with-location-data" out the other end! – TIF Jan 25 at 11:48

I suggest you use gdal rasterize to convert shp to raster. You can use gdal cli or python wrappers, both have very good documentation. Then again use gdal python to read the raster into a numpy array. When you open the raster with gdal python you have access to all of the needed spatial information, store them into variables and use them later when you are writing the simulation results to raster on disk. Numpy arrays themself do not operate with spatial information.

I find gdal and gdal python somewhat tedious to install so if you are not in for a long term processing of spatial data, an easier to install and simpler to use library called rasterio might work for you. I believe it also uses gdal in the background but the installation is simpler.


I've found the easiest way to do this is to install the free, open-source software QGIS, import my map as a layer and use the built-in rasterize tool which worked very quickly. https://www.qgis.org/en/site/

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.