plotting grid files with python and matplotlib

[Ben Buse]

Hi

Here's a simple script to plot a grd file using matplotlib

Matplotlib is nice because can use Philippe Pinard scale bar feature https://github.com/ppinard/matplotlib-scalebar and you can zoom both on the plot and on the color scale bar to change the color levels

I started using gdal which was difficult to get working because I had to pip install the whl for the correct python version from https://www.cgohlke.com/#gdal but then I couldn't get numpy and gdal to let me readasarray stating 'A module that was compiled using NumPy 1.x cannot be run in NumPy 2.0'. So switched to using rasterio. Hence the gdal stuff commented out.

Code Select Expand

#import os
#from osgeo import gdal
import rasterio
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar
import numpy as np
import tkinter
import tkinter.filedialog
#ds = gdal.Open(r"C:\.. .grd")
#geotransform = ds.GetGeoTransform
#rows = ds.RasterYSize
#cols = ds.RasterXSize
imagefileselection=tkinter.filedialog.askopenfile(title="select .grd file to plot",filetypes=(('grd files','grd'),))
src = rasterio.open(imagefileselection.name)
data = src.read()
profile = src.profile

#ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,geotransform[3]+geotransform[5]*rows,geotransform[3]])
#xmingrd = geotransform[0]
#xmaxgrd = geotransform[0]+(geotransform[1]*cols)
#fig.show()
#ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,geotransform[3]+geotransform[5]*rows,geotransform[3]])

grdxpixels = profile['width']
grdypixels = profile['height']

xmingrd = profile['transform'][2]
xmaxgrd = profile['transform'][2]+(profile['transform'][0]*grdxpixels)
ymingrd = profile['transform'][5]+profile['transform'][4]*grdypixels
ymaxgrd = profile['transform'][5]

#transpose changes the axes of 3 dimensional array - as gave error eg Invalid shape (1, 397, 530) for image data - could use try function, and only apply if error
#arr_t = np.transpose(data, (2,1,0))
arr_t = np.transpose(data, (1,2,0))
#plt.show()
fig, ax = plt.subplots()
#ax.imshow(arr_t)
ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,ymingrd,ymaxgrd])
#plt.imshow(arr_t,extent=[xmingrd,xmaxgrd,ymingrd,ymaxgrd])
#scalebar = ScaleBar(profile['transform'][0],"mm",length_fraction=0.5)
scalebar = ScaleBar(1,"mm",length_fraction=0.5)
ax.add_artist(scalebar)
fig.colorbar(ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,ymingrd,ymaxgrd]))
#ax.add_artist(scalebar)
fig.canvas.draw()
print(scalebar.info)
#fig.show()
plt.show()

Here's some screenshots of output

Plot grid file


Zoom on colourbar


Zoom on map

[Ben Buse]

And here's a variant to plot to multiple windows all the oxide.grd files or all the quant.grd files in a folder, to allow each one to be adjusted (in terms of zoom on map or colourbar) before selecting save

Code Select Expand

import os
#from osgeo import gdal
import rasterio
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar
import numpy as np
import tkinter
import tkinter.filedialog
import tkinter.messagebox
answer = tkinter.messagebox.askyesno("Oxide or Elemental","Do you want oxide maps? (If select no elemental maps plotted)")
if answer == True:
    fileend="Oxide.grd"
else:
    fileend="Quant.grd"
folder_dir = tkinter.filedialog.askdirectory(title="Select folder containing .grd files to plot")
#imagefileselection=tkinter.filedialog.askopenfile(title="select .grd file to plot",filetypes=(('grd files','grd'),))
#plotcounter=0
#d={}
for images in os.listdir(folder_dir):
    if (images.endswith(fileend)):
        src = rasterio.open(folder_dir+'/'+images)
        data = src.read()
        profile = src.profile

        #ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,geotransform[3]+geotransform[5]*rows,geotransform[3]])
        #xmingrd = geotransform[0]
        #xmaxgrd = geotransform[0]+(geotransform[1]*cols)
        #fig.show()
        #ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,geotransform[3]+geotransform[5]*rows,geotransform[3]])

        grdxpixels = profile['width']
        grdypixels = profile['height']

        xmingrd = profile['transform'][2]
        xmaxgrd = profile['transform'][2]+(profile['transform'][0]*grdxpixels)
        ymingrd = profile['transform'][5]+profile['transform'][4]*grdypixels
        ymaxgrd = profile['transform'][5]

        #transpose changes the axes of 3 dimensional array - as gave error eg Invalid shape (1, 397, 530) for image data - could use try function, and only apply if error
        #arr_t = np.transpose(data, (2,1,0))
        arr_t = np.transpose(data, (1,2,0))
        #plt.show()
       
        fig, ax = plt.subplots()
        #ax.imshow(arr_t)
        ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,ymingrd,ymaxgrd])
        ax.set_title(images)
        #plt.imshow(arr_t,extent=[xmingrd,xmaxgrd,ymingrd,ymaxgrd])
        #scalebar = ScaleBar(profile['transform'][0],"mm",length_fraction=0.5)
        scalebar = ScaleBar(1,"mm",length_fraction=0.5)
        ax.add_artist(scalebar)
        fig.colorbar(ax.imshow(arr_t,extent=[xmingrd,xmaxgrd,ymingrd,ymaxgrd]))
        #ax.add_artist(scalebar)
        fig.canvas.draw()
        print(scalebar.info)
        #fig.show()
plt.show()