Help getting started with julia geo

[Michael Borregaard]

Hi,
I have migrated my data analysis from R to Julia, with the exception of all the GIS functionality. It appears that Julia has support for shapefiles, a wrapper to GEOS (but I cannot find one to GDAL?) and rudimentary raster in RasterIO, but it is hard to see where to start. Is the geo-functionality in a state where it is usable, and how to get started? Thanks!

Michael

[Yeesian Ng]

It depends on how friendly you expect the packages to be. Julia's support for GIS functionality is nowhere close to that of R's, but is in a good place if you're looking for just reader/writer packages.

There has been a few attempts at wrapping GDAL so far, and I think https://github.com/visr/GDAL.jl's is the most promising. It expects the user to be conversant with the GDAL Data Model, and to work with pointers/handles to them within Julia, and is meant for other packages (e.g. https://github.com/yeesian/GDALUtils.jl) to build upon. These are still fairly alpha-stage packages, so see the tests for example-usage/documentation, and use them at your own risk/etc.

Warnings aside, I think https://github.com/visr/GDAL.jl does not shoehorn users into a single philosophy, and always allows upstream-packages/users the option of bypassing the wrapper/convenience functionality to access the raw C API (which may/may-not be a big deal for you at this point).

The state of support for geographic projections in Julia is also improving, and I'd recommend https://github.com/FugroRoames/Proj4.jl and https://github.com/awbsmith/Geodesy.jl if you're looking for libraries.

[Chris Foster]

On Sat, Apr 30, 2016 at 12:28 AM, Yeesian Ng <ngye...@gmail.com> wrote:
> It depends on how friendly you expect the packages to be. Julia's support
> for GIS functionality is nowhere close to that of R's, but is in a good
> place if you're looking for just reader/writer packages.
>
> There has been a few attempts at wrapping GDAL so far, and I think
> https://github.com/visr/GDAL.jl's is the most promising. It expects the user
> to be conversant with the GDAL Data Model, and to work with pointers/handles
> to them within Julia, and is meant for other packages (e.g.
> https://github.com/yeesian/GDALUtils.jl) to build upon. These are still
> fairly alpha-stage packages, so see the tests for
> example-usage/documentation, and use them at your own risk/etc.
>
> Warnings aside, I think https://github.com/visr/GDAL.jl does not shoehorn
> users into a single philosophy, and always allows upstream-packages/users
> the option of bypassing the wrapper/convenience functionality to access the
> raw C API (which may/may-not be a big deal for you at this point).
>
> The state of support for geographic projections in Julia is also improving,
> and I'd recommend https://github.com/FugroRoames/Proj4.jl and
> https://github.com/awbsmith/Geodesy.jl if you're looking for libraries.

Just a note of warning here - Geodesy is in a state of flux right now,
and after several prototypes we're still not entirely sure what the
high level interface should look like. The main sticking point is
that we'd like to carry CRS information in the point types, but this
can give the API a very intrusive/heavyweight feeling, and in some
circumstances results in an explosion of distinct point types. I've
promised to give a talk about it at juliacon 2016 though, so I guess
that means we'll have to get our act together by then!

The upshot is that the Geodesy API is likely to change a lot in the
near term. It will probably go somewhat in the direction of
awbsmith/Geodesy.jl, but I'm hoping for something a bit more
lightweight. Finding the "right" abstractions here is turning out to
be a rather tricky design problem.

The Proj4.jl package is likely to stay similar to what it is now as
it's mostly just a wrapper around the proj.4 C API. Some of the high
level API will probably change when we figure out how Geodesy.jl and
CoordinateTransformations.jl are going to look.

Cheers,
~Chris

[Michael Krabbe Borregaard]

Thanks a lot, that is really useful. And I guess the answer is it is not ready yet for someone at my user level. I installed the packages and tried running the examples from the GDALUtils readme, but didn't make it past the second line (the package does not appear to be exporting any methods).

I need it for things like 1. load two shapefiles, get them in the same projection, make a new shape based on their intersection and plot it in colors by some attribute; 2. load a raster, use the shapefile to mask/crop it, plot it, and calculate the sum and standard deviation of values.

I appreciate all your hard work on this, and will look greatly forward to the day I can let totally go of my umbilical cord to R :-)

--
You received this message because you are subscribed to the Google Groups "julia-geo" group.
To unsubscribe from this group and stop receiving emails from it, send an email to julia-geo+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Yeesian Ng]

Yup, I don't think the libraries are stable/usable enough for end-users in their current state either. Maybe/optimistically in 1 to 2 years.

I need it for things like 1. load two shapefiles, get them in the same projection, make a new shape based on their intersection and plot it in colors by some attribute; 2. load a raster, use the shapefile to mask/crop it, plot it, and calculate the sum and standard deviation of values.

I'll be grateful if you can help provide some worked examples (in R/etc). It will be really helpful for developers to know what tasks to target, and improve the API for common/idiomatic approaches.

[Michael Krabbe Borregaard]

Sure here is a worked example in R - is this useful for you?


library(rgdal) #readOGR, spTransform, CRS, proj4string
library(rgeos) #gIntersection, gArea

#Example 1: Get the land area in equal area grid cells

# calculates the intersection between two shapefiles
get.area <- function(gr, cst){ #gr and wm are just object isd
  int <- gIntersection(cst, grid[grid$WorldMapID==gr,],byid=TRUE)
  if (length(int)==0) return(0)
  gArea(int)
}

# read a shapefile with a 1x1 deg equal area grid (e.g. NCEAS EASE)
grid <- readOGR("../Grid shapefile", "EqualAreaGrid")
wm.ids <- grid@data$WorldMapID

# read the world coastline
coast <- readOGR("../Grid shapefile", "World_coast")
# put on the same projection
coast <- spTransform(coast,CRS(proj4string(grid)))

# Area of landmass in each grid cell in this projection
landareas <- sapply(wm.ids, function(gr) get.area(gr, coast))



# Example 2 - extract climatic variables in a given area
# for convenience I use each of the land mass polygons from before

library(raster) #raster, crop, mask, extent, getValues

#load the raster and set the NA value
NPP <- raster("../Climatic data/MOD17A3_Science_NPP_mean_00_14.tif")
NAvalue(NPP) <- 65535

#Extract sub rasters for each polygon
landmass_list <- lapply(1:nrow(coast), function(cst){
  tmp <- crop(rast, extent(coast[cst,]))
  tmp <- mask(tmp, coast[cst,])
  tmp
})

# Calculate the mean NPP for all landmasses
meanNPP <- sapply(landmass_list, function(rst) mean(getValues(rst), na.rm = T))
names(meanNPP) <- coast$Name

# plot the polygons with a color based on NPP
cols <- terrain.colors(64)
zlim <- range(meanNPP, na.rm = T)
colorvec <- meanNPP - zlim[1]
colorvec <- floor(colorvec * (length(cols) - 1)/(zlim[2] - zlim[1])) +
  1
plot(coast, col = cols[colorvec])

--