Traduire EPSG: 4326 coordonnées de latitude lon en xy en utilisant pyproj

Aug 19 2020

J'essaie de travailler sur les détails de l'utilisation de pyproj pour traduire les coordonnées lat lon du flux USGS Earthquake en coordonnées xy correspondant à ces positions sur un fond de carte d'une largeur et d'une hauteur déterminées et d'une origine en haut à gauche. Jusqu'à présent, mon code utilise naïvement PIL / oreiller pour traduire les positions lon lat en xy et dessiner des cercles sans prendre en compte la projection EPSG: 4326 du fond de carte:

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)

J'ai parcouru la documentation et j'ai le sentiment général que j'utiliserais quelque chose comme:

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

avec des dimensions de largeur et de hauteur. L'exécution de la ligne ci-dessus donne l'erreur:

TypeError: p1 must be a Proj class

Mise à jour J'obtiens des valeurs pour:

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

Dans ce cas, une entrée (lon, lat) de (-110, 39)donne la sortie 1.923, 0.688qui n'est pas correcte. Ce dont j'ai besoin, c'est d'une fonction de projection qui prend le lon, la latitude et la largeur et la hauteur de l'image et génère les coordonnées x, y. J'écrirais cela moi-même, mais je ne trouve pas d'équation appropriée pour cette conversion et j'espérais trouver cette capacité dans proj.

MISE À JOUR Voici le fond de carte que j'utilise.

Et voici le dessin qui en résulte. Je suppose que les points centraux des cercles marquant les tremblements de terre sont désactivés, car j'utilise des coordonnées cartésiennes pour lon, lat, mais le fond de carte a été créé avec une projection et doit donc avoir une distorsion à partir des coordonnées cartésiennes x, y.

De plus, je réalise que je dois mettre ma légende sur une échelle logarithmique. Je n'y suis pas encore arrivé.

Réponses

5 FSimardGIS Aug 22 2020 at 04:12

Lorsque je superpose votre fond de carte sur une projection équirectangulaire, il s'adapte parfaitement. Cela signifie donc que votre fond de carte utilise simplement la simple "projection" de tracer les longitudes comme x et les latitudes comme y. Ce qui signifie que la seule chose dont vous devez vous soucier est de vérifier la taille de l'image, ce que vous faites déjà dans vos formules.

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

Ceux-ci fonctionneront correctement pour vos besoins, au pixel près. En fait, j'ai vérifié les tremblements de terre sur l'USGS à cette date, et vos points sont tous au bon endroit!