Add point coordinates to a GeoTable

Specific Domains Geo
,
1

Hello everyone,

I am trying to learn Julia and move some of my geostats projects to this language. I found that the GeoStats package is well regarded, hence I am using it. I have a grid and I want to separate the east from the west side of it. However, I am struggling to find how to retrieve x coordinates and append these values in the geotable. It follows a reproducible example:

using CSV, DataFrames, GeoStats

Aᵢ = CSV.read("a.csv", DataFrame; header = false) # 900x1 dataframe
grid = CartesianGrid(30, 30)
Aᵢ = georef(Aᵢ, grid)
f(x) = ustrip(coords(centroid(x)).x)<15 

Aᵢ = @transform(Aᵢ, :geometry => ByRow(f) => :Country) # ERROR: LoadError: ArgumentError: invalid expression syntax

Importantly, the f function is working as intended. This can be seen if I do

f.(Aᵢ.geometry)

Any help would be highly appreciated!

1 Like
2

Hi @phchavesmaia , welcome to the forum!

Try the geosplit utility:

  geosplit(object, fraction, [normal])

  Split geospatial object into two parts where the first part has a fraction of
  the elements. Optionally, the split is performed perpendicular to a normal
  direction.

I believe you want

geosplit(geotable, 0.5, (1, 0))
1 Like
3

Thank you very much, @juliohm , this indeed solved part of my problem!

Now, I just wanted to create a column (flag) in the main GeoTable indicating which cells belong to each section. What I am trying is:

using CSV, DataFrames, GeoStats, CairoMakie

data = CSV.read("a.csv", DataFrame; header = false)
rename!(data,[:Aᵢ])
data.cell = 1:size(data,1) 
grid = CartesianGrid(30, 30)
data = georef(data, grid)
data.Country = data.cell.∈ Ref(geosplit(data,0.5,(1,0))[1].cell) # ERROR: type GeoTable has no field Country

Any advice?

4

It feels like an instance of the XY problem. Could you please describe what are trying to achieve as your final objective? Did you read the GDSJL book?

https://juliaearth.github.io/geospatial-data-science-with-julia

5

One possible way to create a new column based on the x coordinate:

@transform(geotable, :column = coords(:geometry).x < 100km)
6

My final objective is to replicate Redding and Rossi-Hansberg (2017) in julia, while the original code is in Matlab. I know that @SebKrantz has already dabbled into this endeavor, but I want to take my own shot at it.

In this paper, the authors generate a 30 by 30 grid representing two bordering countries that trade with each other while facing export tarifs. I am trying to generate the such grid, indicating to which country a cell belongs, in order to reproduce their findinds.

That said, differently from them, I would like to reproduce the model while using georeferenced data, which is the kind of data I would actually use in a real research scenario.

7

I tried to do it, but it returns the following error: ArgumentError: column name “x” not found in the data frame; existing most similar names are: “Aᵢ”

What I did was to create an auxiliary table to store the values (as in the GDSJL book). I made all the operations regarding values using this auxiliary table and then just plugged the geometries back. Somewhat convoluted, but it worked!

8

I believe we encountered a bug in which the expression coords(:geometry).x is not correctly understood by the @transform macro. Will take a look, thanks for reporting.

Please let us know if something else is not working as expected, or if you need additional help.

1 Like
9

Here is a valid solution:

fun(point) = coords(point).x < 10km

@transform(geotable, :column = fun(:geometry))

Try to define the function outside the macro.

1 Like
10

@phchavesmaia we identified the root of the problem: the expression obj.field is converted to obj.:field by Julia and the resulting symbol :field is incorrectly interpreted as a column in the @transform macro.

We will treat this as a special case like other packages do (e.g., DataFramesMeta.jl). You can watch the corresponding issue to get notifications when it is fixed:

1 Like
11

This worked perfectly, thank you!

1 Like