Skip to content

Instantly share code, notes, and snippets.

@hagerty
Created September 14, 2016 18:01
Show Gist options
  • Save hagerty/68f622c79486d746cfc042e42befa84b to your computer and use it in GitHub Desktop.
Save hagerty/68f622c79486d746cfc042e42befa84b to your computer and use it in GitHub Desktop.
Distance transform from geoJSON to a bumpy array
import json
from shapely.geometry import Polygon, shape, Point
import gdal
import numpy as np
import sys
fn = sys.argv[1]
path = sys.argv[2]
def Pixel2World ( geoMatrix, i , j ):
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
return(1.0 * i * xDist + ulX, -1.0 * j * xDist + ulY)
ds8 = gdal.Open(path+'8band/'+'8band_'+fn+'.tif')
ds3 = gdal.Open(path+'3band/'+'3band_'+fn+'.tif')
geoTrans = ds8.GetGeoTransform()
with open(path + 'vectorData/geoJson/'+fn+'_Geo.geojson','r') as f:
js = json.load(f)
dist = np.zeros((ds8.RasterXSize, ds8.RasterYSize))
for i in range(ds8.RasterXSize):
for j in range(ds8.RasterYSize):
point = Point(Pixel2World( geoTrans, i , j ))
pd = -100000.0
for feature in js['features']:
polygon = shape(feature['geometry'])
newpd = point.distance(polygon.boundary)
if False == polygon.contains(point):
newpd = -1.0 * newpd
if newpd > pd :
pd = newpd
dist[i,j] = pd
np.save(path+'CosmiQ_distance/'+fn+'.distance',dist)
@Imaget
Copy link

Imaget commented Mar 2, 2017

Great and thank you very much for providing code. Can you please help me how to show dist ?

@bquast
Copy link

bquast commented Jun 26, 2017

As it is now, it contains 5 indentation errors, fixed version here:

https://gist.github.com/bquast/eec2b10fe2594c166e2bfd2c5d67dde5

@StefankinaOlya
Copy link

Hello. Can somebody show me a json file?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment