Methods to extract pixel coordinates after plotting

Hello,

I plot rotated geological formation polygons outputted by another tool called pyGPlates (coordinates in longitude and latitude). I would like to have a way to know where the plotted polygons are located in terms of pixel coordinates (relative to the generated image). I wonder if there is a way to know such information while using Fig.plot() or any other similar functions.

Example:

Above is the image I plotted using Fig.plot(). I would like to know the pixel coordinates of those colorful polygons in the middle.

Thanks for any help in advance.

Hi @xdhqyz!! Welcome to the GMT Forum.

Do you have a file of that colorful polygons? Could you share here?

Hello Esteban,

Thanks for the reply. Here’s a file that was generated by pyGPlates. The system does not allow uploading of .gmt file format so I renamed it to .txt. Please change it back to .gmt when running. reconstructed_geom.gmt.txt (21.8 KB)

I think that you could try something like this to get the X and Y positions (in cm) of each point of the polygon.

gmt mapproject -Rg reconstructed_geom.gmt.txt -JW15c

Just change the value of R and J to your values.

To get the pixels instead, I think that you would need to multiply by the resolution of your image (Pixels per cm).
Be careful with the frame and anotations outside the map.

Thank you. Since we are using PyGMT, may I ask if pygmt.project is the equivalent function of gmt mapproject?

No, that is the equivalent of gmt project.

So is there no equivalent of mapproject in pygmt?

The pyGMT users/developers can answer that - I cannot.

Okay I see. Thank you anyways

Toe extend @Esteban82’s example to give integer pixels, you would need something like this (again -R -J on you) and of course do it in pyGMT:

echo 10 10 | gmt mapproject -Rg -JW15c | gmt math STDIN -Ca 300 MUL RINT =

This says that lon = 10, lat = 10 has pixel coordinates

145 1278

I have a doubt. Why 300? Is not 300 the resolution in pixels per inch? Should we use the resolution of pixel per centimeters instead?

Wel, this is all tied into PROJ_LENGTH_UNIT. If inch then 300 if dpi is 300, else if cm then dry/2.54.

mapproject is not wrapped in PyGMT yet. There is already a feature request (https://github.com/GenericMappingTools/pygmt/issues/1986) but we haven’t found time to implement it. In other words, if you like, you can help us implement it.