Question about Distance limit

[ForVerm]

I start to work with geosardine. Could someone tell me what are the units for distance limit and power? I assumed the distance limit was in the same unit as the resolution but it doesn't seems to be right. Thanks for your help.

[sahitono]

Hi, @ForVerm could you please kindly show how you executed, the result, and the expected results?

[ForVerm]

Hi Sahitono, thanks for your help.
I am not a coder, only a final "tester". Somebody else put your code in form for our utilisation. I manually test few distance limit values as well as power values. My resolution is 15 arcsec express in degrees (0.00416667). Extent is in Lat/Long. I have a point cloud surface which is my input file. Points are 50 to 1500 meters appart. I am trying to get a result as close as possible to this reference file. My better fit is with distance limit = 200 and power = 4. However I don't understand what this 200 stands for. Are they meters? Power is a kind of smoothing, righ?

Parameters are in a json file.
{"Table_SEP": "St-Laurent_Golfe_2019.csar",
"coin_inf_lat": 46.81250090,
"coin_inf_lon": -71.22083710,
"coin_sup_lat": 52.28805500,
"coin_sup_lon": -55.77066520,
"resolution_lat": 0.00416667,
"resolution_lon": 0.00416667,
"distance_limit": 200,
"power": 4,
"epsg": 8237,
"longlat_distance": "harvesine",
"Output": "StLauGolf_15arcsec_d200p4.txt"}

.py follow (descriptions in french)

import coverage_utils as cu
import numpy as np
import geosardine as dine
import rasterio
import fiona
import json
import os
import sys

#On spécifie un path
APPLICATION_PATH = ''
if getattr(sys, 'frozen', False):
APPLICATION_PATH = os.path.dirname(sys.executable)
elif file:
APPLICATION_PATH = os.path.dirname(file)

#On importe les paramètres inscrits dans le fichier de configuration de type JSON
with open(os.path.join(APPLICATION_PATH, 'config_IDW.json'), 'r') as fich:
param = fich.read()
params = json.loads(param)

#On initialise la resolution spatiale en degrées décimales puisque qu'on fournit de coordonnées géodésique

#(voir documentation Geosardine epsg=4326)

reso_lat = params['resolution_lat']
reso_lon = params['resolution_lon']

#On initialise l'extent de la grille qu'on veut créer avec des coordonnées Lat et Lon
#(j'ai choisi le coin inférieur gauche et le coin supérieur droit)
coin_inf_lon = params['coin_inf_lon']
coin_inf_lat = params['coin_inf_lat'] - reso_lat #on ajuste la latitude pour inclure le premier point de l'extent dans l'interpolation
coin_sup_lon = params['coin_sup_lon']
coin_sup_lat = params['coin_sup_lat']

#On initialise la distance limite de à être interpolé (distance_limit : float, default=0, maximum distance to be interpolated, can't be negative)
distance_limit = params['distance_limit']

#Paramètre en relation au rayon d'influence des points adjascents (power : float, default=2, how smooth the interpolation will be)
power = params['power']

epsg = params['epsg']

method used to calculate distance in spherical / ellipsoidal

* If harvesine, calculation will be faster but less accurate

* If vincenty, calculation will be slower but more accurate

longlat_distance = params['longlat_distance']

#On importe le module coverage_utils
table_sep = cu.CoverageUtils()

#On crée un objet nuage contenant notre table de séparation en format point clouds
nuage = table_sep.get_coverage(os.path.join(APPLICATION_PATH, params['Table_SEP']))

#On itère sur notre nuage de points pour aller chercher les coordonnées de chacuns des points
block = table_sep.iterate_cloud_blocks(nuage, 'XYZ')

#On initialise les variable XY et Z qui serviront d'intrant à Geosardine-IDW (les coord. XY et Z)
XY = block[0][0:2]
Z = block[0][2]

#On fait une boucle pour mettre nos coordnnées sous le format nd.array

qui convient en input à Geosardine-IDW

for i in block[1::]:
XY = np.vstack((XY, i[0:2]))
Z = np.append(Z, i[2])

#Vérifications que les dimensions concordent
print('Nombre de points importées du fichier CSAR')
print(Z.shape)
print(XY.shape)

#On fait l'interpollation avec Geosardine-IDW (Voir documentation--> https://sahitono.github.io/geosardine/index.html)
interpolated = dine.interpolate.idw(XY, Z, spatial_res=(reso_lat, reso_lon), epsg=epsg,
longlat_distance=longlat_distance,
extent=(coin_inf_lon, coin_inf_lat, coin_sup_lon, coin_sup_lat),
power=power,
distance_limit=distance_limit)

#On trouve le nombre de colonnes et de lignes qu'il y a dans l'extent qu'on a spécifié en fonction de la résolution
count_lon = (coin_sup_lon - coin_inf_lon) // reso_lon
count_lat = (coin_sup_lat - coin_inf_lat) // reso_lat

#Vérification
print('Dimensions de la grille en longitude')
print(count_lon)
print('Dimension de le grille en latitude')
print(count_lat)
print('Nombre de valeurs interpollées')
print(interpolated.array.size)

#On itère en fonction de la résolution pour créer notre output
XYZ_out = []

#Après vérificarion on constate que interpolated.array est donnée lignes par lignes en descendant à partir du coin en haute à gauche jusqu'en bas
#Donc il faut créer manuellement notre gille XY en fonction de l'extent et de la résolution. Par ailleurs, portons attention à l'énoncé suivant
#-->(int(count_lat)-i) qui décale d'une latitude la grille vers le haut pour faire concorder l'interpollation IDW.
for i in range(int(count_lat-1), -1, -1):
for j in range(int(count_lon)):
XYZ_out.append((coin_inf_lat + (int(count_lat)-i)reso_lat, coin_inf_lon + jreso_lon, interpolated.array[i][j]))

#On fait le trie sur les latitudes en ordre décroissant afin de fournir les coordonnées XYZ dans le bon format pour

l'utilitaire de conversion vers .bin Geoidconverter (https://eye4software.com/hydromagic/)

sorted_par_latitudes = sorted(XYZ_out, key=lambda tup: tup[0], reverse=True)

#On écrit dans un fichier texte les coord XYZ lequel on pourra convertir avec l'utilitaire GeoideCOnverter
#(https://eye4software.com/hydromagic/)
with open(params['Output'], 'w') as fp:
fp.write('\n'.join('{} {} {}'.format(format(round(x[1], 8), '.8f'),
format(round(x[0], 8), '.8f'),
format(round(x[2], 3), '.3f'))
for x in sorted_par_latitudes))

print('Votre fichier ' + params['Output'] + ' est créé')

[ForVerm]

Sorry for this big bold text. I really don't know where it comes from ???

[sahitono]

that is correct , distance_limit is the distance limit that has the same unit as the reference system and power is the smoothness degree.

it looks like you are using epsg 8327 which is a projected reference system. Therefore, your coordinate and resolution should be in meter. epsg is only used for defining the source coordinate.
so, if your coordinate is latlong WGS84, you should use epsg 4326.

[ForVerm]

Thanks for your help!

[sahitono]

you're welcome. you can try the latest version as it supports all epsg for non projected epsg.

I'll mark it as answered