length(Bkde.x) must be equal to size(Bkde.density, 2) or size(Bkde.density, 2) + 1. Same for y. Difference is grid registration, a must know feature that is not GMT specific but valid for all grids. See 4. Standardized command line options — GMT 6.5.0 documentation
Like this?
G = mat2grid(Float32.(Bkde.density), x=(Bkde.density, 2), y=(Bkde.density, 2), proj4=“+proj=moll +lon_0=0 +x_0=0 +y_0=0”)
Changing the resolution to 1km with Rasters.jl is a single command…
rast_1k = resample(rast; res=1000, crs=your_crs_in_meters)
I’m not sure why you would use GMT.jl for that when you already had a Raster.
(But yes just going directly to 1000 is better than resampling - But you would just wrap the output of KernelDensity as a Raster?)
My fault, but after a couple of days and the problem not solved I decided to offer an alternative.
Yeah, I’ve never used KernelDensity.jl so can’t really help with that part.
But writing a raster is the easy bit, probably just one line after you work out getting a matrix from the Interpolations.jl pdf.
Like:
pdf_matix = ... # whatever you need to do with KernelDensity.jl / Interpolations.jl
write("some.tif", Raster(pdf_matrix, (X(xs), Y(ys)); crs))
We could probably add a KernelDensity.jl extension to Rasters.jl so these things are easier, and it just takes any Table/DataFrame and outputs a kernel density Raster.
Hi @Raf
Thank you
I eventually ended up with rasterizing the output of the pdf function (see code below)
I still have two questions
Looking at the line you just wrote you specifiy Raster(pdf_matrix, (X(xs), Y(ys)) instead of Y(reverse(ys)) …
I also did not reverse the y coordinate and I do not know why this was suggested by @Fliks can anyone explain me?
Also, after masking the raster and writing it into the disk I noticed that the missing values are transformed into -Inf
This is the warning message I got
Warning: missing cant be written with gdal, missinval for Float64 of -Inf used instead
I woud like to substitute them with NA but I can pass only strings to replace_missing
So ideally I should replace them with “NA”
The thing is that I would like to import them in R afterwards where missing values are coded with NA …
what do you suggest?
Thank you very much for all suppport
# Extract x and y columns
x_values = df.x
y_values = df.y
# Extract x and y coordinates
coordinates = hcat(x_values, y_values)
bw=30000::Int64
# Defining the bounding of the World ESRI:54009 Mollweide
# and resolution of 1kmx1km (1000m x 1000m)
res=1000
xmin=-17601618
xmax=17601617
ymin=-9018991
ymax=8751339
lengthx=round(Int,(xmax + abs(xmin)) / res)
lengthy=round(Int,(ymax + abs(ymin)) / res)
rx = KernelDensity.kde_range((xmin,xmax), lengthx)
ry = KernelDensity.kde_range((ymin, ymax), lengthy)
# Create a bivariate kernel density estimator
Bkde = kde(coordinates,(rx,ry), bandwidth=(bw,bw))
# Create a raster file
ras=Raster(Bkde.density, (X(rx), Y((ry))), crs=EPSG(54009))
# Upload the World Shapefile
World=Shapefile.Table("World.shp") |> DataFrame
plot(World.geometry)
# Crop the raster using the World Shapefile
ras=Rasters.crop(ras;to=World.geometry)
# Mask the raster using the World Shapefile
ras_masked=Rasters.mask(ras,with=World.geometry)
Plots.plot(ras_masked)
# Define the coords reference system
crs = ArchGDAL.toWKT(ArchGDAL.importPROJ4("+proj=moll +lon_0=0 +x_0=0 +y_0=0"))
ras_masked=setcrs(ras_masked,crs)
mean(skipmissing(ras_masked))
Rasters.write("SpatialKernelall.tif", ras_masked)
The standard in gdal/most of the spatial world is for Y to be reversed.
Rasters used to do this for you automatically, but it turned out to be annoying in some edge cases of skew and rotation and it means round-trip read/write returned a different array. So now you need to reverse it yourself if you want it. e.g. with reverse(raster, Y) afterwards, or just make the Y axis backwards from the start.
For missing values, Julia doesn’t have the concept of NA and GDAL/tif doesn’t either. It has either sentinel values (the most common) or a separate bitmask layer for in some formats.
Rasters only uses a sentinel value: missingval which you can set manually in replace_missing or just let it choose typemin(eltype(rast)). But R should read that fine and convert them to Na, its part of the generic gdal interface that R packages also use.
replace_missing(rast, newmissingval) should work? I don’t understand what ou mean with:
I would like to substitute them with NA but I can pass only strings to replace_missing
There is no NA in julia. Do you mean you want to use NaN ? That is only possible with Float64.
Also, you don’t need to do this conversion, Rasters/GDAL will do it automatically for you in write:
crs = ArchGDAL.toWKT(ArchGDAL.importPROJ4("+proj=moll +lon_0=0 +x_0=0 +y_0=0"))
ras_masked=setcrs(ras_masked,crs)
You can just do crs = ProjString("+proj=moll +lon_0=0 +x_0=0 +y_0=0").
Hi @Raf you are right and R automatically consider them as Na
Thank you so much for your support
PS I agree, an extension of KernelDensity to Rasters may make things easier 
Since you have learned how it works now, feel free to add it in a PR!
Basically it will just be your code in a kerneldensity(data; to, res, size, crs, missingval) function like most other Rasters methods. data could be any Table.jl compatible table or GeoInterface.jl compatible PointTrait or MultiPointTrait or vector of them.
We would add it in an extension based on KernelDensity.jl so we don’t add the dependencies for everyone. There are other extensions in Rasters.jl you can copy, (like for gdal, with e.g. resample).
If you get the bones of it in place and working I can help get it finished.