Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import shapefile
- import matplotlib.pyplot as plt
- import math
- import pycrs
- import pyproj
- import sys
- def toCartesian(p):
- degToRad=math.pi/180
- return (math.cos(p[1]*degToRad)*math.cos(p[0]*degToRad),math.cos(p[1]*degToRad)*math.sin(p[0]*degToRad),math.sin(p[1]*degToRad))
- def fromCartesian(p):
- degToRad=math.pi/180
- return (math.atan2(p[1],p[0])/degToRad,math.asin(p[2])/degToRad)
- def normalize(p):
- dis=math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2])
- return (p[0]/dis,p[1]/dis,p[2]/dis)
- def cross(a,b):
- return (a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0])
- def dist(a,b):
- lat1=a[1]
- lat2=b[1]
- lon1=a[0]
- lon2=b[0]
- p=math.pi/180
- a=0.5 - math.cos((lat2-lat1)*p)/2 + math.cos(lat1*p) * math.cos(lat2*p) * (1-math.cos((lon2-lon1)*p))/2 # https://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula
- return 2*math.asin(math.sqrt(a))
- def nearPointGreatCircle(a,b,c): # https://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
- a1,b1,c1=map(toCartesian,(a,b,c))
- G=cross(a1,b1)
- F=cross(c1,G)
- t=normalize(cross(G,F))
- return fromCartesian(t)
- def onSegment(a,b,c):
- return abs(dist(a,b)-dist(a,c)-dist(c,b))<10**-9
- def segDist(a,b,c):
- t=nearPointGreatCircle(a,b,c)
- if onSegment(a,b,t):
- return dist(c,t)
- return min(dist(c,a),dist(c,b))
- koordGrad=open('koordList.txt','rb').read().split(b'\n')
- stanov={}
- zupGrada={}
- popis=open('popis2011.txt','rb').read().split(b'\n') # uzimamo koji je grad u kojoj županiji iz popisa stanovništva; prvo smo tražili težinsku udaljenost gdje je težina svakog grada/općine bio njezin broj stanovnika, ali budući da je bilo čudno da se udaljenost od mora mijenja kroz vrijeme na kraju nismo koristili broj stanovnika
- for x in popis:
- z=x.split(b'\t')
- if z[1]==b'':
- continue
- if z[4][:5]!=b'Vinod':
- imeGrada=z[1]+b' '+z[4]
- else:
- imeGrada=z[4]
- if ' – '.encode() in imeGrada:
- imeGrada=imeGrada.split(' – '.encode())[0]
- stan=int(b''.join(z[-1].split(b',')))
- stanov[imeGrada]=stan
- zupGrada[imeGrada]=z[0]
- for x in koordGrad:
- z=x.split(b'\t')
- imeGrada=z[2].split(b'\r')[0]
- xCoord=float(z[1])
- yCoord=float(z[0])
- zup=zupGrada[imeGrada]
- pop=stanov[imeGrada]
- sf=shapefile.Reader('eea_v_3035_100_k_coastline-poly_1995-2017_p_v03_r00.zip')
- fromcrs=pycrs.load.from_file("eea_v_3035_100_k_coastline-poly_1995-2017_p_v03_r00/EEA_Coastline_20170228.prj")
- fromproj=pyproj.crs.CRS.from_wkt(fromcrs.to_ogc_wkt())
- tocrs=pycrs.parse.from_epsg_code(4326)
- toproj=pyproj.crs.CRS.from_wkt(tocrs.to_ogc_wkt())
- transformer=pyproj.Transformer.from_crs(fromproj,toproj)
- print(sf)
- fig,ax=plt.subplots()
- hrv_obala=sf.shape(33540).points[328000:350000] # dio obale koji je relevantan za hrvatsku
- minx=min([hrv_obala[i][0] for i in range(len(hrv_obala))])
- maxx=max([hrv_obala[i][0] for i in range(len(hrv_obala))])
- miny=min([hrv_obala[i][1] for i in range(len(hrv_obala))])
- maxy=max([hrv_obala[i][1] for i in range(len(hrv_obala))])
- for i in range(len(hrv_obala)):
- hrv_obala[i]=transformer.transform(hrv_obala[i][0],hrv_obala[i][1]) # pretvorimo projekciju koja se koristila u koordinate
- xs=[p[0] for p in hrv_obala]
- ys=[p[1] for p in hrv_obala]
- ax.plot(xs,ys)
- otoci=[]
- for i in range(70972):
- boundingBox=sf.shape(i).bbox
- if boundingBox[0]>minx and boundingBox[1]>miny and boundingBox[2]<maxx and boundingBox[3]<maxy: # ako je otok potpuno sadržan u bounding boxu hrvatske obale onda je to hrvatski otok kojeg također treba gledati
- otok=[]
- for x in sf.shape(i).points:
- otok.append(transformer.transform(x[0],x[1]))
- xs=[p[0] for p in otok]
- ys=[p[1] for p in otok]
- ax.plot(xs,ys)
- otoci.append(otok)
- plt.show()
- totStanovnika={}
- zbrUdalj={}
- zupanije=[]
- for x in koordGrad:
- z=x.split(b'\t')
- imeGrada=z[2].split(b'\r')[0]
- xCoord=float(z[1])
- yCoord=float(z[0])
- zup=zupGrada[imeGrada]
- pop=stanov[imeGrada]
- pop=1
- toc=(xCoord,yCoord)
- minDist=-1
- for i in range(len(hrv_obala)-1):
- di=segDist(hrv_obala[i],hrv_obala[i+1],toc)
- if minDist==-1 or di<minDist:
- minDist=di
- for obala_otoka in otoci:
- for i in range(len(obala_otoka)-1):
- di=segDist(obala_otoka[i],obala_otoka[i+1],toc)
- if minDist==-1 or di<minDist:
- minDist=di
- minDist*=6371 # radijus zemlje u km
- print(imeGrada,xCoord,yCoord,zup,pop,minDist)
- if zup not in zupanije:
- zupanije.append(zup)
- totStanovnika[zup]=0
- zbrUdalj[zup]=0
- totStanovnika[zup]+=pop
- zbrUdalj[zup]+=pop*minDist
- ff=open('res.txt','wb')
- for x in zupanije:
- prosUdalj=zbrUdalj[x]/totStanovnika[x]
- print(x,prosUdalj,totStanovnika[x])
- ff.write(x+b': '+str(prosUdalj).encode()+b'\n')
- ff.close()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
