4 min to read
Satellite Images - easy to work with
I have been in the business of analysing satellite images for over 2 years now. I and my collaborators code primarily in Python for our back-end. There are several Python libraries which make our lives much easier to manupulate geographical data.
Here is the list of these libraries:
- Rasterio - useful for plotting
- Fiona - useful for coordinate system discovery
- Shapefile - useful for loading satellite image data
- Geopandas - builds upon pandas, shapely, descartes and matplotlib - very useful for data manipulation similar to pandas & coordinate system manipulaton.
Now let’s showcase what they can do in an example scenario.
Imagine that we have downloaded 1 tif (tif stands for geotiff) image which contains the NVDI indices calculated based on a particular satellite image. NDVI stands for Normalized Difference Vegetation Index and it is the go-to vegetation index for developers and geospatial specialists when inspecting the state of vegetation on Earth. The index is calculated by using the values of the red and near-infra red radiation captured by the satellite image.
A Typical NDVI image looks something like this:
The darker the green colour is, the denser the vegetation is. NVDI takes values between 0 and 1. A commonly used rule of thumb says that when the NVDI value is smaller than 0.2 - what we see on the image is bare soil or soil with a some grass. When the values are smaller than 0.4, we usually look at bushes, and when the values are greather than 0.5 - we see a forest (-.8 for a very dense forest).
Find out what the coordinate system of the image is
When a tif image is opened in Python, one of the first things we want to know is what the coordinate system of the satellite image is. Fiona comes in handy in here:
import fiona import shapefile as shp path_to_image = 'path/image.tif' c = fiona.open(path_to_image) original_crs = c.crs['init'] print(original_crs) >> 'epsg:3763'
Convert an image to a different coordinate system
epsg:3763 is one of the coordinate systems available. A Standard coordinate system is WSG84 (WSG stands for (World Geodetic System](https://en.wikipedia.org/wiki/World_Geodetic_System)) and 84 is its laest revision. The
espg code for WSG84 is
espg:4326. We can see that the example image I chose to load is based on a different coordinate system so it is good to convert it to WSG84. used geopandas to do this:
import geopandas as gpd gpd_shape = gpd.read_file(path_to_image) gpd_shape.crs = original_crs['init'] # we set the CRS we have identified when we loaded the image using fiona gpd_shape = gpd_shape.to_crs("EPSG:4326")
Export the image to GeoJSON format
path_to_image_geojson= 'path/image.geojson' gpd_shape.to_file(path_to_image_geojson, driver='GeoJSON')
Find the Bouding Box of shapes
_____ separate this blog post into 2 - 1 for shapes, and 2 for satellite images, and 3- putting it all together
bbox = shape.bbox >> [-47903.63338843786, -265251.9380933815, -31310.261804120528, -257247.52285356744]
Write Bounding Box to a file
upper_left = [bbox, bbox] upper_right = [bbox, bbox] lower_left = [bbox, bbox] lower_right = [bbox, bbox] # write bbox to a shape file (just to visualize in Google Earth) w = shp.Writer(path_to_aglomerados_bbox) # assign POLYGON (5) as shapeType w.shapeType = shp.POLYGON # add a "name" field of type "Character" w.field('BBOX', 'C') w.poly([[upper_left, # upper left upper_right, # upper right lower_right, # lower right lower_left, # lower left upper_left]]) # first point again w.record('Shape_Bbox') w.close()
FYI - I used this Linux terminal command to create the GIF of two PNG images I had generated:
convert -delay 60 -loop 0 img1.png img2.png result.gif
Calculate and Classify Normalized Difference Results with EarthPy
Handling the error TypeError: unhashable type: ‘MultiPolygon’
Expanding MultiPolygons to Polygons in GeoPandas