将 scipy 插值图转换为 .tiff 文件并保存到目录

Question

我已经使用 scipy.interpolate 模块生成了插值图。我需要一些帮助将地图保存为 .tiff 文件并将其保存到我的目录中。但是，我不确定是否需要将它转换为 numpy 数组，因为它需要在每个单元格中具有纬度、经度和插值数据。任何帮助将不胜感激！

这是数据。 nutrition.csv 文件 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

使用此代码生成的地图需要保存为地理参考 .tiff 文件。

Answer 1

这是我之前回答过的问题的后续。要将数组保存到 geotiff，您需要确定地理变换，这意味着您需要知道数组左上角的坐标以及 x 和 y 的分辨率。

对于您的数据，它可能是这样工作的：

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)

我没有使用 holoviews 所以我无法测试它，您可能必须先更改 numpy 数组中的变量或使用不同的代码来定义 width、dtype 等.根据您的数据集中北的定义方式，转换可能会有所不同，例如

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

还要检查 xres 和 yres 的符号。

查看 rasterio rasterio.readthedocs.io 的文档 - 它非常简洁明了，将帮助您充分理解，使这项工作适合您。

显然，如果您的数据在不同的投影中，您必须以不同的方式定义 crs，但我相信 latlong 是您所需要的。

将 scipy 插值图转换为 .tiff 文件并保存到目录

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

python

numpy

scipy

pandas