Я создал карту интерполяции с помощью модуля scipy.interpolate
. Мне нужна помощь, чтобы сохранить карту в виде файла .tiff
и сохранить его в моем каталоге. Однако я не уверен, нужно ли мне преобразовывать его в массив numpy или нет, так как в каждой ячейке должны быть широта, долгота и интерполированные данные. Любая помощь приветствуется!
Вот данные. nutrition.csv
файл можно найти здесь.
#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
.
Итак, это продолжение вашего вопроса, на который я ответил ранее. Чтобы сохранить массив в геотиф, вам нужно определить геотрансформацию, а значит, вам нужно знать координаты левого верхнего угла вашего массива и разрешение по 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 — это то, что вам нужно.
Отлично, большое спасибо за вашу помощь! Это очень помогает.