[GPlates-discuss] Reconstruct rasters with pyGPlates

[Simon Williams]

Hi Ben,

pygplates does not currently have an official raster capability that matches the same capability in GPlates, in the sense that you cannot simply specify a GPlates-readable raster file, load it in to a feature collection, do stuff with it with the same speed as GPlates using the GPU.

Also, for the purposes of this reply I’ll assume that when you say ‘deformed’, you mean reconstructed based on rigid plate motions rather than a reconstruction model with actual ‘deforming’ regions. If you do mean the latter, then can stop reading now…..

Instead, there does exist some experimental code built on top of pygplates to make raster reconstruction more efficient than it would be by simply treating each grid node as a point (or multipoint).

[ One of the main issues for this kind of calculation is that the step where each grid node is assigned a plateid would be much slower in python compared to how it is done within GPlates. The atom_utils repository contains some functions that John wrote to speed up this process using spatial trees, which typically makes a big difference. ]

John may have a better (more well documented) example, the best one I can think of can be found here:
https://protect-au.mimecast.com/s/1CRcCgZowLHBMZrAiNAYi5?domain=github.com

See specifically the function ‘reconstruct_raster’

An example usage would be:

    Raster = GplatesRaster(filename)

    (reconstructed_point_lons,
     reconstructed_point_lats,
     reconstructed_zvalues) = reconstruct_raster(raster, static_polygons, rotation_model,
                                             from_time, to_time,
                                             grid_sampling=2.0)

    plt.scatter(reconstructed_point_lons,unreconstructed_point_lats,

                c=reconstructed_zvalues)

Where static_polygons is a list of polygon features required for the plateid assignments, rotation_model is an associated set of rotations to make these polygons move, from_time would be zero for a raster in present-day coordinates and to_time would be the desired reconstruction time.

Note that:
- This code is not exactly working on 2D numpy arrays, rather a feature of the (even more experimental) GPlates raster class, which is assumed to be coming from a netcdf file, see here https://protect-au.mimecast.com/s/JgGPCjZrzqHOWQ2nCRrM7t?domain=github.com. The class contains arrays of the lats and longs of the grid nodes as well as an array of the z values. You could modify the first few lines of the reconstruct_raster function to avoid using this class.
- The code returns arrays of reconstructed longs and lats together with z values, but not on a regular lat-long grid - you’d need to to resample as an extra step.
- if you try to use this code and hit unfound dependency problems then feel free to let me know.

Cheers,
Simon








On 16 Dec 2018, at 11:22 am, Ben Mather <ben.mather at sydney.edu.au<mailto:ben.mather at sydney.edu.au>> wrote:

Hi team!

I was wondering if it is possible to reconstruct raster images back in time like it is in regular GPlates? It would be very useful to, say, have a 2D numpy array that could be deformed with respect to a rotation model.

Cheers,
Ben

Dr. Ben Mather | Computational Geophysicist
Room 418, Madsen Building F09
School of Geoscience, Faculty of Science
The University of Sydney, NSW 2006
m +61 422 470 117
w benmather.info<https://protect-au.mimecast.com/s/ppjXCk8vAZtwxz7OHQboEj?domain=benmather.info>
t twitter.com/BenRMather<https://protect-au.mimecast.com/s/x3BIClxwB5CMQNy2fq9_6T?domain=twitter.com>


_______________________________________________
GPlates-discuss mailing list
GPlates-discuss at mailman.sydney.edu.au<mailto:GPlates-discuss at mailman.sydney.edu.au>
https://mailman.sydney.edu.au/mailman/listinfo/gplates-discuss

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.sydney.edu.au/pipermail/gplates-discuss/attachments/20181217/cc15eb24/attachment-0001.html>