Перевод EPSG: 4326 долготных координат в xy с помощью pyproj

Aug 19 2020

Я пытаюсь проработать детали использования pyproj для перевода координат широты и долготы из канала USGS Earthquake в координаты xy, соответствующие этим положениям на базовой карте определенной ширины и высоты и начала координат в верхнем левом углу. Мой код до сих пор наивно использует PIL / подушку для перевода долготы в координаты xy и рисования кругов без учета проекции EPSG: 4326 базовой карты:

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)

Я просмотрел документацию и понял, что могу использовать что-то вроде:

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

с некоторыми размерами ширины и высоты. Выполнение строки выше дает ошибку:

TypeError: p1 must be a Proj class

Обновление. Я получаю значения для:

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

В этом случае ввод (lon, lat) (-110, 39)дает 1.923, 0.688неправильный вывод . Мне нужна функция проекции, которая принимает длину, широту, ширину и высоту изображения и выводит координаты x, y. Я бы написал это сам, но я не могу найти подходящего уравнения для этого преобразования и надеялся найти эту возможность в proj.

ОБНОВЛЕНИЕ Вот базовая карта, которую я использую.

А вот и получившийся рисунок. Я предполагаю, что центральные точки для кругов, отмечающих землетрясения, отключены, потому что я использую декартовы координаты для lon, lat, но базовая карта была создана с помощью проекции и, следовательно, должна иметь искажение от декартовых координат x, y.

Кроме того, я понимаю, что мне нужно поместить легенду в логарифмическую шкалу. Еще не дошел до этого.

Ответы

5 FSimardGIS Aug 22 2020 at 04:12

Когда я накладываю вашу базовую карту на равнопрямоугольную проекцию, она идеально подходит. Это означает, что ваша базовая карта на самом деле просто использует простую «проекцию» отображения долготы как x и широты как y. Это означает, что единственное, о чем вам нужно беспокоиться, это проверить размер изображения, что вы уже делаете в своих формулах.

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

Они будут работать правильно для ваших нужд с точностью до ближайшего пикселя. Фактически, я проверил землетрясения на USGS в тот день, и все ваши очки находятся в правильном месте!