Extract grid value at x, y location

CompulsoryCoffee · October 25, 2022, 3:22pm

Hi,

I tried many things without success and I am keen to stop installing more libraries

I loaded a grid

G = gmtread("/GIS/g2015_pga.tif")
G
100×100×1 Raster{Float32,3} with dimensions: 
  X Projected LinRange{Float64}(-1.7754e5, 6.10362e5, 100) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected LinRange{Float64}(3.37598e6, 2.91423e6, 100) ReverseOrdered Regular Intervals crs: WellKnownText,
  Band Categorical 1:1 ForwardOrdered
with missingval: 0.0

I have a CSV file with columns for x coordinates and y coordinates (± 7 millions Float points)

How can I extract the Raster value for all those points?

Thanks

Joaquim · October 25, 2022, 4:15pm

CompulsoryCoffee:

G = gmtread("/GIS/g2015_pga.tif")
G
100×100×1 Raster{Float32,3} with dimensions: 
  X Projected LinRange{Float64}(-1.7754e5, 6.10362e5, 100) ForwardOrdered Regular Intervals crs: WellKnownText,

This G was not read with GMT.gmtread was it? Because that’s not how it pretty prints data. If I open a Geotiff grid, it prints (an exampe)

julia> G = gmtread("fract.tif")
title: Grid imported via GDAL
Pixel node registration used
x_min: -12.009784735812133      x_max :-1.9902152641878672      x_inc :0.019569471624266144     n_columns :512
y_min: 34.99021526418787        y_max :45.00978473581213        y_inc :0.019569471624266144     n_rows :512
z_min: -14.40057373046875       z_max :27.675704956054688
PROJ: +proj=longlat +datum=WGS84 +no_defs
512×512 Matrix{Float32}:
 -0.857803    -0.956818   -1.04156   -1.0193     …   0.51535     0.0821075  -0.322647   -0.562424   -0.751602

What does typeof(G) show in your case?

Anyway, if you have a GMTgrid type (my G) and your CSV file is simple (in the sense that it has only x,y,z rows) you’d do

D1 = gmtread("fille.csv");
D2 = grdtrack(G, D);

CompulsoryCoffee · October 26, 2022, 2:39pm

Thanks Joaquim (again)

Sorry I didn’t copy the full print from the raster - it seems OK:

julia> G
title: Grid imported via GDAL
Pixel node registration used
x_min: -177509.3592     x_max :618390.6407999999        x_inc :100.0    n_columns :7959
y_min: 2.9141275377e6   y_max :3.3886275377e6   y_inc :100.0    n_rows :4745
z_min: 0.0      z_max :0.7368712425231934
PROJ: +proj=utm +zone=45 +datum=WGS84 +units=m +no_defs
4745×7959 Matrix{Float32}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

I have an issue with the CSV file I think

x,y

322100.813851,3081391.061107

322108.793984,3081388.668287

158210.720118,3063161.547234

95122.282451,3049126.951143

95122.293548,3049138.132305

When I tried to load it

julia> pts
BoundingBox: [-177472.06767, 617812.376218, 2.914211052193e6, 3.36171579452e6]
Error showing value of type GMTdataset{Float64, 2}:
ERROR: MethodError: Cannot `convert` an object of type String to an object of type Symbol

Not sure what is happening there (I tried without the header but it didn’t work)…

More than often my CSV have a lot of attributes attached to them that I load in a DataFrame (df).

Can I do :

pts = Tuple(hcat(df.x, df.y))
values = grdtrack(pts, G)

???

Joaquim · October 26, 2022, 3:38pm

GMT.jl version 0.43.1 just got merged so please update first. The error

ERROR: MethodError: Cannot `convert` an object of type String to an object of type Symbol

is coming from prettytables and you can ignore it, but it should be gone in latest version. If I put those 5 points in a file (with or without a #x,y comment line) it works fine (ignore also the warning. Must see why it’s happening)

julia> gmtread("aa.dat")
┌ Warning: Failed to parse file aa.dat for file_has_time!(). Error was: BoundsError(SubString{String}["322100.813851,3081391.061107"], (2,))
└ @ GMT C:\Users\joaqu\.julia\dev\GMT\src\gmtreadwrite.jl:229
Comment:        ["x,y"]
BoundingBox: [95122.282451, 322108.793984, 3.049126951143e6, 3.081391061107e6]
5×2 GMTdataset{Float64, 2}
 Row │         col.1      col.2
     │       Float64    Float64
─────┼──────────────────────────
   1 │     3.22101e5  3.08139e6
   2 │     3.22109e5  3.08139e6
   3 │     1.58211e5  3.06316e6
   4 │ 95122.3        3.04913e6
   5 │ 95122.3        3.04914e6

For more complex CSV files you may need to use the CSV/DataFrames packages but at the end you must pass a simple Matrix (not a Tuple) to the grdtrack module. It accepts both matrices or GMTdataset types.

CompulsoryCoffee · October 28, 2022, 11:28am

Everthing was working fine but I am having issues with a bunch of rasters and I am not sure why (EPSG is 32645, it might cause an issue):

**julia>** G = gmtread("su30m_nepal_propa_17843_nb_propagations.tif")

**ERROR:** UndefVarError: data not defined

Stacktrace:

[1] **get_image(**API::Ptr{Nothing}, object::Ptr{Nothing}**)**

@ GMT ~/.julia/packages/GMT/mzT4h/src/gmt_main.jl:631

[2] **GMTJL_Get_Object(**API::Ptr{Nothing}, X::GMT.GMT_RESOURCE**)**

@ GMT ~/.julia/packages/GMT/mzT4h/src/gmt_main.jl:883

[3] **gmt(**::String**)**

@ GMT ~/.julia/packages/GMT/mzT4h/src/gmt_main.jl:326

[4] **gmtread(**fname::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}**)**

@ GMT ~/.julia/packages/GMT/mzT4h/src/gmtreadwrite.jl:148

[5] **gmtread(**fname::String**)**

@ GMT ~/.julia/packages/GMT/mzT4h/src/gmtreadwrite.jl:68

[6] top-level scope

@ REPL[15]:1

edit:

I had this information on Windows as well:

Warning 1: PROJ: proj_create_from_database: C:\Program Files\AppJ\PostgreSQL\13\share\contrib\postgis-3.1\proj\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata. It comes from another PROJ installation.
Warning 1: The definition of projected CRS EPSG:32645 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
Warning 1: PROJ: proj_create_from_database: C:\Program Files\AppJ\PostgreSQL\13\share\contrib\postgis-3.1\proj\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata. It comes from another PROJ installation.
Warning 1: The definition of projected CRS EPSG:32645 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

Joaquim · October 28, 2022, 1:18pm

It looks like it’s trying to open your file as an image instead of grid (we don’t use much the word raster here because it mixes these two different types)

Two things:

Can you make that file available so I can try to debug the failures cause?
try G = gdlaread( "su30m_nepal_propa_17843_nb_propagations.tif")

CompulsoryCoffee · October 28, 2022, 2:27pm

Here it is:

Joaquim · October 28, 2022, 5:16pm

So your file is BIGGGGGGGGGGG

julia> gdalread("su30m_nepal_propa_26516_nb_propagations.tif")
Gridline node registration used
x_min: -177634.86776953447      x_max :618425.1322304655        x_inc :10.0     n_columns :46661
y_min: 2.914137533405694e6      y_max :3.380737533405694e6      y_inc :10.0     n_rows :79607
z_min: 0.0      z_max :31671.0
PROJ: +proj=utm +zone=45 +datum=WGS84 +units=m +no_defs
79607×46661 Matrix{UInt32}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

79607 * 46661 * 4 / 1024/1024/1024 ~= 13.8 GB

As you see I could still read it with gdalread but being of that size it looks not very usable (and it actually more or less froze my computer).

Will later try to find time why gmtread fails to read it but my immediate guess is that it’a a consequence of the data type (UInt32 (??)).

CompulsoryCoffee · November 7, 2022, 9:12am

Yes it is. I received the format from an external party so I wouldn’t know much about the processing unfortunately.

Could you gdalread with a coordinates box as constrain?

Thanks for giving it a shot.

Joaquim · November 7, 2022, 11:50am

help?> gdalread
search: gdalread gdalrasterize gdalwrite gdaltranslate gdalvectortranslate

  gdalread(fname::AbstractString, opts=String[]; gdataset=false, kwargs...)

  Read a raster or a vector file from a disk file and return the result either as a GMT type (the default) or a GDAL
  dataset.

    •  fname: Input data. It can be a file name, a GMTgrid or GMTimage object or a GDAL dataset

    •  opts: List of options. The accepted options are the ones of the gdal_translate utility. This list can be in
       the form of a vector of strings, or joined in a simgle string.

    •  gdataset: If set to true forces the return of a GDAL dataset instead of a GMT type.

    •  kwargs: This options accept the GMT region (-R) and increment (-I)

  Returns
  –––––––––

  A GMT grid/image or a GDAL dataset

So a limits=(....) option is supposed to work (though I have some distant recalls of issues the increment option)