Convert scipy interpolation map to .tiff file and save to directory-CodePudding

I have generated an interpolation map using the scipy.interpolate module. I am needing some help saving the map as a .tiff file and saving it to my directory. However, I'm not sure if I need to convert it to a numpy array or not, as it needs to have the latitude, longitude, and the interpolated data in each cell. Any help would be much appreciated!

Here is the data. The nutrition.csv file can be found here.

#Import modules
import pandas as pd
import numpy as np
import os
import shapely
import geopandas as geo
import glob
import holoviews as hv
hv.extension('bokeh')
from scipy.interpolate import griddata, interp2d
import fiona
import gdal
import ogr

#Read file
nut = pd.read_csv('nutrition.csv') #Data to be interpolated

#Minimum and maximum longtude values
lon_min = nut['longitude'].min()
lon_max = nut['longitude'].max()

#Create range of nitrogen values
lon_vec = np.linspace(lon_min, lon_max, 30) #Set grid size

#Find min and max latitude values
lat_min = nut['latitude'].min()
lat_max = nut['latitude'].max()

# Create a range of N values spanning the range from min to max latitude
# Inverse the order of lat_min and lat_max to get the grid correctly
lat_vec = np.linspace(lat_max,lat_min,30,)

# Generate the grid
lon_grid, lat_grid = np.meshgrid(lon_vec,lat_vec)

#Cubic interpolation
points = (nut['longitude'], nut['latitude'])
targets = (lon_grid, lat_grid)
grid_cubic = griddata(points, nut['n_ppm'], targets, method='cubic', fill_value=np.nan)

#Generate the graph
map_bounds=(lon_min,lat_min,lon_max,lat_max)


map_cubic = hv.Image(grid_cubic, bounds=map_bounds).opts(aspect='equal',
                                                         colorbar=True,
                                                         title='Cubic',
                                                         xlabel='Longitude',
                                                         ylabel='Latitude',
                                                         cmap='Reds')

map_cubic

The map produced with this code needs to be saved as a georeferenced .tiff file.

CodePudding user response：

So this is the follow up of your question that I answered earlier. To save an array to a geotiff you need to determine the geotransform, which means you need to know the coordinates of the upper left corner of your array and the resolution in x and y.

For your data it might work like this:

xres = lon_vec[1]-lon_vec[0]
yres = lat_vec[1]-lat_vec[0]

from rasterio.transform import Affine

transform = Affine.translation(lon_vec[0] - xres / 2, lat_vec[0] - yres / 2) * Affine.scale(xres, yres)

with rasterio.open(
    '/tmp/new.tif',
    'w',
    driver = 'GTiff',
    height = map_cubic.shape[0],
    width = map_cubic.shape[1],
    count = 1,
    dtype = map_cubic.dtype,
    crs = ' proj=latlong',
    transform = transform,
) as dst:
    dst.write(map_cubic, 1)

I'm not using holoviews so I can't test it, you might have to change your variable in a numpy array first or use different code to define the width, dtypeetc. Depending on how north is defined in your data set the transform might be different, e.g.

transform = Affine.translation(lon_vec[0] - xres / 2, lat_vec[-1] - yres / 2) * Affine.scale(xres, yres)

Also check the sign of xres and yres.

Check out the documentation for rasterio rasterio.readthedocs.io - it's really nice and concise and will help you understand enough to make this work for you.

Obviously, you have to define crs differently if your data is in a different projection, but I believe latlong is what you need.