Skip to content

Commit 3f352fb

Browse files
gabrielaznargabrielaznar
committed
read only one channel + refactor
1 parent cd11d06 commit 3f352fb

File tree

2 files changed

+185
-38
lines changed

2 files changed

+185
-38
lines changed

‎climada/entity/exposures/nightlight.py‎

Lines changed: 68 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -245,55 +245,85 @@ def load_nightlight_nasa(bounds, req_files, year):
245245
# 0, 1, 2, 3 longitude range for min and max longitude
246246
in_lon_nb = math.floor(in_lon[0]/21600), math.floor(in_lon[1]/21600)
247247

248-
prev_idx = -1
248+
nightlight = sparse.lil.lil_matrix([])
249+
idx_info = [0, -1, False] # idx, prev_idx and row added flag
249250
for idx, file in enumerate(BM_FILENAMES):
251+
idx_info[0] = idx
250252
if not req_files[idx]:
251253
continue
252254

253-
aux_nl = Image.open(path.join(SYSTEM_DIR, file.replace('*', str(year))))
254-
aux_nl = sparse.csc.csc_matrix(aux_nl.getchannel(0))
255-
# flip X axis
256-
aux_nl.indices = -aux_nl.indices + aux_nl.shape[0] - 1
257-
aux_nl = aux_nl.tolil()
255+
with Image.open(path.join(SYSTEM_DIR, file.replace('*', str(year)))) \
256+
as im_nl:
257+
aux_nl = im_nl.getchannel(0)
258+
cut_nl_nasa(aux_nl, idx_info, nightlight, in_lat, in_lon,
259+
in_lat_nb, in_lon_nb)
258260

259-
aux_bnd = []
260-
if idx/2 % 4 == in_lon_nb[0]:
261-
aux_bnd.append(int(in_lon[0] - (idx/2%4)*21600))
262-
else:
263-
aux_bnd.append(0)
261+
idx_info[1] = idx
264262

265-
if idx % 2 == in_lat_nb[0]:
266-
aux_bnd.append(in_lat[0] - ((idx+1)%2)*21600)
267-
else:
268-
aux_bnd.append(0)
263+
coord_nl[0, 0] = coord_nl[0, 0] + in_lat[0]*coord_nl[0, 1]
264+
coord_nl[1, 0] = coord_nl[1, 0] + in_lon[0]*coord_nl[1, 1]
269265

270-
if idx/2 % 4 == in_lon_nb[1]:
271-
aux_bnd.append(int(in_lon[1] - (idx/2%4)*21600) + 1)
272-
else:
273-
aux_bnd.append(21600)
266+
return nightlight.tocsr(), coord_nl
274267

275-
if idx % 2 == in_lat_nb[1]:
276-
aux_bnd.append(in_lat[1] - ((idx+1)%2)*21600 + 1)
277-
else:
278-
aux_bnd.append(21600)
279-
280-
if prev_idx == -1:
281-
nightlight = sparse.lil.lil_matrix((aux_bnd[3]-aux_bnd[1],
282-
aux_bnd[2]-aux_bnd[0]))
283-
nightlight = aux_nl[aux_bnd[1]:aux_bnd[3], aux_bnd[0]:aux_bnd[2]]
284-
elif idx%2 == prev_idx%2:
285-
nightlight = sparse.hstack([nightlight, \
286-
aux_nl[aux_bnd[1]:aux_bnd[3], aux_bnd[0]:aux_bnd[2]]])
287-
else:
288-
nightlight = sparse.vstack([nightlight, \
289-
aux_nl[aux_bnd[1]:aux_bnd[3], aux_bnd[0]:aux_bnd[2]]])
268+
def cut_nl_nasa(aux_nl, idx_info, nightlight, in_lat, in_lon, in_lat_nb,
269+
in_lon_nb):
270+
"""Cut nasa's nightlight image piece (1-8) to bounds and append to final
271+
matrix.
290272
291-
prev_idx = idx
273+
Parameters:
274+
aux_nl (PIL.Image): nasa's nightlight part (1-8)
275+
idx_info (list): idx (0-7), prev_idx (0-7) and row_added flag (bool).
276+
nightlight (sprse.lil_matrix): matrix with nightlight that is expanded
277+
in_lat (tuple): min and max latitude indexes in the whole nasa's image
278+
in_lon (tuple): min and max longitude indexes in the whole nasa's image
279+
in_lat_nb (tuple): for min and max latitude, range where they belong
280+
to: upper (0) or lower (1) row of nasa's images.
281+
on_lon_nb (tuple): for min and max longitude, range where they belong
282+
to: 0, 1, 2 or 3 column of nasa's images.
283+
"""
284+
idx, prev_idx, row_added = idx_info
285+
aux_nl = sparse.csc.csc_matrix(aux_nl)
286+
# flip X axis
287+
aux_nl.indices = -aux_nl.indices + aux_nl.shape[0] - 1
292288

293-
coord_nl[0, 0] = coord_nl[0, 0] + in_lat[0]*coord_nl[0, 1]
294-
coord_nl[1, 0] = coord_nl[1, 0] + in_lon[0]*coord_nl[1, 1]
289+
aux_bnd = []
290+
if int(idx/2) % 4 == in_lon_nb[0]:
291+
aux_bnd.append(int(in_lon[0] - (int(idx/2)%4)*21600))
292+
else:
293+
aux_bnd.append(0)
295294

296-
return nightlight.tocsr(), coord_nl
295+
if idx % 2 == in_lat_nb[0]:
296+
aux_bnd.append(in_lat[0] - ((idx+1)%2)*21600)
297+
else:
298+
aux_bnd.append(0)
299+
300+
if int(idx/2) % 4 == in_lon_nb[1]:
301+
aux_bnd.append(int(in_lon[1] - (int(idx/2)%4)*21600) + 1)
302+
else:
303+
aux_bnd.append(21600)
304+
305+
if idx % 2 == in_lat_nb[1]:
306+
aux_bnd.append(in_lat[1] - ((idx+1)%2)*21600 + 1)
307+
else:
308+
aux_bnd.append(21600)
309+
310+
if prev_idx == -1:
311+
nightlight.resize((aux_bnd[3]-aux_bnd[1], aux_bnd[2]-aux_bnd[0]))
312+
nightlight[:, :] = aux_nl[aux_bnd[1]:aux_bnd[3], aux_bnd[0]:aux_bnd[2]]
313+
elif idx%2 == prev_idx%2 or prev_idx%2 == 1:
314+
# append horizontally in first rows e.g 0->2 or 1->2
315+
nightlight.resize((nightlight.shape[0],
316+
nightlight.shape[1] + aux_bnd[2]-aux_bnd[0]))
317+
nightlight[:aux_bnd[3]-aux_bnd[1], -aux_bnd[2]+aux_bnd[0]:] = \
318+
aux_nl[aux_bnd[1]:aux_bnd[3], aux_bnd[0]:aux_bnd[2]]
319+
else:
320+
# append vertically in lasts rows and columns e.g 0->1 or 2->3
321+
if not row_added:
322+
nightlight.resize((nightlight.shape[0] + aux_bnd[3] - aux_bnd[1],
323+
nightlight.shape[1]))
324+
idx_info[2] = True
325+
nightlight[-aux_bnd[3]+aux_bnd[1]:, -aux_bnd[2]+aux_bnd[0]:] = \
326+
aux_nl[aux_bnd[1]:aux_bnd[3], aux_bnd[0]:aux_bnd[2]]
297327

298328
def unzip_tif_to_py(file_gz):
299329
""" Unzip image file, read it, flip the x axis, save values as pickle

‎climada/entity/exposures/test/test_nighlight.py‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@
33
"""
44
import unittest
55
import numpy as np
6+
from scipy import sparse
7+
68
from climada.entity.exposures import nightlight
79
from climada.util.constants import SYSTEM_DIR
810

@@ -71,6 +73,121 @@ def test_download_nightlight_files(self):
7173

7274
# The same length but not the correct length
7375
self.assertRaises(ValueError, nightlight.download_nl_files, (1, 0, 1), (1, 1, 1))
76+
77+
def test_cut_nl_nasa_1_pass(self):
78+
"""Test cut_nl_nasa situation 2->3->4->5."""
79+
nl_mat = sparse.lil.lil_matrix([])
80+
in_lat = (21599, 21600)
81+
in_lon = (43199, 43200)
82+
# 0 2 4 6 (lat: Upper=0) (lon: 0, 1, 2, 3)
83+
# 1 3 5 7 (lat: Lower=1) (lon: 0, 1, 2, 3)
84+
in_lat_nb = (1, 0)
85+
in_lon_nb = (1, 2)
86+
87+
idx_info = [2, -1, False]
88+
aux_nl = np.zeros((21600, 21600))
89+
aux_nl[21599, 21599] = 100
90+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
91+
in_lon, in_lat_nb, in_lon_nb)
92+
93+
self.assertEqual(nl_mat.shape, (1, 1))
94+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
95+
96+
idx_info[0] = 3
97+
idx_info[1] = 2
98+
aux_nl = np.zeros((21600, 21600))
99+
aux_nl[0, 21599] = 101
100+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
101+
in_lon, in_lat_nb, in_lon_nb)
102+
103+
self.assertEqual(nl_mat.shape, (2, 1))
104+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
105+
self.assertEqual(nl_mat.tocsr()[1, 0], 101.0)
106+
107+
idx_info[0] = 4
108+
idx_info[1] = 3
109+
aux_nl = np.zeros((21600, 21600))
110+
aux_nl[21599, 0] = 102
111+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
112+
in_lon, in_lat_nb, in_lon_nb)
113+
114+
self.assertEqual(nl_mat.shape, (2, 2))
115+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
116+
self.assertEqual(nl_mat.tocsr()[1, 0], 101.0)
117+
self.assertEqual(nl_mat.tocsr()[0, 1], 102.0)
118+
119+
idx_info[0] = 5
120+
idx_info[1] = 4
121+
aux_nl = np.zeros((21600, 21600))
122+
aux_nl[0, 0] = 103
123+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
124+
in_lon, in_lat_nb, in_lon_nb)
125+
126+
self.assertEqual(nl_mat.shape, (2, 2))
127+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
128+
self.assertEqual(nl_mat.tocsr()[1, 0], 101.0)
129+
self.assertEqual(nl_mat.tocsr()[0, 1], 102.0)
130+
self.assertEqual(nl_mat.tocsr()[1, 1], 103.0)
131+
132+
def test_cut_nl_nasa_2_pass(self):
133+
"""Test cut_nl_nasa situation 3->5."""
134+
nl_mat = sparse.lil.lil_matrix([])
135+
in_lat = (21599, 21599)
136+
in_lon = (43199, 43200)
137+
# 0 2 4 6 (lat: Upper=0) (lon: 0, 1, 2, 3)
138+
# 1 3 5 7 (lat: Lower=1) (lon: 0, 1, 2, 3)
139+
in_lat_nb = (1, 1)
140+
in_lon_nb = (1, 2)
141+
142+
idx_info = [3, -1, False]
143+
aux_nl = np.zeros((21600, 21600))
144+
aux_nl[0, 21599] = 100
145+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
146+
in_lon, in_lat_nb, in_lon_nb)
147+
148+
self.assertEqual(nl_mat.shape, (1, 1))
149+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
150+
151+
idx_info[0] = 5
152+
idx_info[1] = 3
153+
aux_nl = np.zeros((21600, 21600))
154+
aux_nl[0, 0] = 101
155+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
156+
in_lon, in_lat_nb, in_lon_nb)
157+
158+
self.assertEqual(nl_mat.shape, (1, 2))
159+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
160+
self.assertEqual(nl_mat.tocsr()[0, 1], 101.0)
161+
162+
def test_cut_nl_nasa_3_pass(self):
163+
"""Test cut_nl_nasa situation 2->4."""
164+
nl_mat = sparse.lil.lil_matrix([])
165+
in_lat = (21600, 21600)
166+
in_lon = (43199, 43200)
167+
# 0 2 4 6 (lat: Upper=0) (lon: 0, 1, 2, 3)
168+
# 1 3 5 7 (lat: Lower=1) (lon: 0, 1, 2, 3)
169+
in_lat_nb = (0, 0)
170+
in_lon_nb = (1, 2)
171+
172+
idx_info = [2, -1, False]
173+
aux_nl = np.zeros((21600, 21600))
174+
aux_nl[21599, 21599] = 100
175+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
176+
in_lon, in_lat_nb, in_lon_nb)
177+
178+
self.assertEqual(nl_mat.shape, (1, 1))
179+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
180+
181+
idx_info[0] = 4
182+
idx_info[1] = 2
183+
aux_nl = np.zeros((21600, 21600))
184+
aux_nl[21599, 0] = 101
185+
nightlight.cut_nl_nasa(aux_nl, idx_info, nl_mat, in_lat,
186+
in_lon, in_lat_nb, in_lon_nb)
187+
188+
self.assertEqual(nl_mat.shape, (1, 2))
189+
self.assertEqual(nl_mat.tocsr()[0, 0], 100.0)
190+
self.assertEqual(nl_mat.tocsr()[0, 1], 101.0)
74191

75192
# Execute Tests
76193
TESTS = unittest.TestLoader().loadTestsFromTestCase(TestNightLight)

0 commit comments

Comments
 (0)