2

I'm trying to work out the details of using pyproj to translate the lat lon coordinates from the USGS Earthquake feed to xy coordinates corresponding to those positions on a basemap of a determined width and height and origin at the upper left. My code so far naively uses PIL/pillow to translate lon lat to xy positions and draw circles without taking into account the EPSG:4326 projection of the basemap:

from PIL import Image
from PIL import ImageDraw

# open basemap image file
basemap = Image.open(basemap_path).convert('RGBA)

# resize to desired map size
basemap.thumbnail(width, height, Image.LANCZOS)

# get proportional height
width_bmp, height_bmp = basemap.size

# create background frame and paste basemap on it
img=Image.new('RGB',(width, height), color = '#000000')
img.paste(basemap, (0,0), basemap)
draw = ImageDraw.Draw(img, 'RGBA')

width_scale = width/360
height_scale = height_bmp/180

# usgs data has been parsed into a list
for quake in earthquake_list:
    lon = float(quake["longitude"])
    lat = float(quake["latitude"])
    mag = float(quake["mag"])

    # want to use pyproj to translate coordinates here instead of the following
    cx = (lon + 180) * width_scale
    cy = (90-lat) * height_scale
    r  = scaleRadius(mag)

    # draw earthquake circles
    draw.ellipse((cx-r, cy-r, cx+r, cy+r), fill = colormap(mag))

draw = ImageDraw.Draw(img)
img.save(filepath, quality=100)

I've gone through the documentation and have a general sense that I would use something like:

cx, cy = pyproj.transform("EPSG:4326", "xy", lon, lat)

with some width and height dimensions. Running the line above gives the error:

TypeError: p1 must be a Proj class

Update I do get values for:

p = Proj(proj = 'longlat', ellps='WGS84')
cx,cy = p(lon, lat)

In this case a (lon, lat) input of (-110, 39) gives the output 1.923, 0.688 which isn't right. What I need is a projection function that takes the lon, lat and the width and height of the image and outputs the x, y coordinates. I would write this on my own, but I can't find a suitable equation for that conversion and was hoping to find this capability in proj.

UPDATE Here is the basemap I am using.

enter image description here

And here is the resulting drawing. I'm assuming the centerpoints for the circles marking earthquakes are off, because I am using Cartesian coordinates for lon, lat, but the basemap was created with a projection and thus must have a distortion from x, y Cartesian coordinates.

enter image description here

Also, I realize I need to put my legend on a logarithmic scale. Haven't gotten to that yet.

9
  • Note, EPSG:4326 coordinates are in lat/long order – nmtoken Aug 19 '20 at 8:45
  • I'm not sure how you can have a (not-georeferenced) global EPSG:4326 basemap. Either the base map is in EPSG:4326 and therefore georeferenced, or it's just a map image with no CRS. Can you edit the question to expand on what you mean by the statement. – nmtoken Aug 19 '20 at 9:05
  • I've rephrased that description. Thanks. – interwebjill Aug 19 '20 at 9:18
  • X/Y in relation to what? – bugmenot123 Aug 20 '20 at 20:11
  • 1
    I have updated the post with images of the basemap as well as the result of drawing the circles at the lon, lat coordinates from the USGS data. – interwebjill Aug 21 '20 at 7:49
5
+50

When I overlay your basemap on an Equirectangular projection, it fits perfectly. So this means your basemap really just uses the simple "projection" of plotting longitudes as x and latitudes as y. Which means the only thing you have to worry about is to check the size of the image, which you already do in your formulas.

cx = (lon + 180) * width_scale
cy = (90-lat) * height_scale

These will work properly for your needs, to the nearest pixel. In fact, I verified the Earthquakes on USGS on that date, and your points are all in the right location!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.