2D Array with a different value for every country [message #94219] |
Wed, 01 March 2017 07:11  |
c.betancourt
Messages: 5 Registered: March 2017
|
Junior Member |
|
|
I would like to fold plumes of air pollutants with emissions of these pollutants (“fold” = hadamard product in this case). My plumes are 2D Arrays, with longitudes and latitudes as coordinates. Since I have emission factors for specific countries, I would like to make another 2D array with the known emission factors. For example all dots within France would have the value 0.55.
Do you have an idea how to do that?
Thank you,
Clara
|
|
|
|
|
|
|
Re: 2D Array with a different value for every country [message #94231 is a reply to message #94230] |
Thu, 02 March 2017 01:27   |
Helder Marchetto
Messages: 520 Registered: November 2011
|
Senior Member |
|
|
On Thursday, March 2, 2017 at 9:56:00 AM UTC+1, c.beta...@fz-juelich.de wrote:
>> ;plumes of air pollutants matrix of 720x1440 elements
>> pap = fltarr(720,1440)
>> ;the scale is 0.25x0.25. Take this into account when plotting.
>> ;now fill the pap array...
>> ...
>> ;now generate the emission array:
>> emiss = fltarr(720,1440)
> -> -> -> emis = ???
>> ;must have the same dimensions and scale of the pap array.
>> ;now fill the emiss array with the emission values
>> ...
>> ;the idea is that all the pixels (indices) having coordinates in France will have a value of 0.55 and so on.
>> ;now do the hadamard product of two
>> result = pap * emiss
>
>
>
>
> Sorry to confuse you. My only problem is that I don´t have an emission array yet. I need something like a world map as an array. A similar problem would be: I would like to plot a world map where france is green, italy is yellow and germany is red. Hope this is precise enough...?
This is not my field, but here you might start looking for an answer:
http://www.harrisgeospatial.com/docs/mappingcontinents.html
The second example seems pretty close to what you want:
; Define a map of Europe.
map = MAP('STEREOGRAPHIC', FILL_COLOR = 'Light Blue', $
LIMIT = [30.0, -15.0, 68.0, 55.0])
; Add the country outlines and fill color.
mc = MAPCONTINENTS(/COUNTRIES, FILL_COLOR='beige')
; Add the rivers.
rivers = MAPCONTINENTS(/RIVERS, COLOR='blue')
Cheers,
Helder
|
|
|
|
Re: 2D Array with a different value for every country [message #94233 is a reply to message #94232] |
Thu, 02 March 2017 03:05   |
Helder Marchetto
Messages: 520 Registered: November 2011
|
Senior Member |
|
|
On Thursday, March 2, 2017 at 11:01:43 AM UTC+1, c.beta...@fz-juelich.de wrote:
> Dear Helder,
>
> OK. Still, there remain two problems: I can choose the color for water and land areas. I did not find any option or keyword to access a certain country. The second problem is, this is not an array that I could modify or fold with anything.
>
> Thank you,
> Clara
Hi,
has I said, countries or imagery of the earth surface and stuff is not my stuff.
But I have some good news: here you can find info about how to get coordinates of territories out of shapefile:
http://www.idlcoyote.com/code_tips/extractpoly.php
By the way, when you draw the maps of Europe, IDL gets the shapefile (where the borders are defined) from this file:
IDL path\resource\maps\shape\cntry08.shp
Now the bad news. The pro cgExtractShape is not for free:
http://www.idlcoyote.com/idldoc/forsale/index.html
You should ask David for a price, might be not that bad... But I sincerely have no idea.
I bought some time ago activecontour.pro and I'm very happy with it.
Cheers,
Helder
|
|
|
Re: 2D Array with a different value for every country [message #94234 is a reply to message #94233] |
Thu, 02 March 2017 03:38   |
Helder Marchetto
Messages: 520 Registered: November 2011
|
Senior Member |
|
|
On Thursday, March 2, 2017 at 12:05:37 PM UTC+1, Helder wrote:
> On Thursday, March 2, 2017 at 11:01:43 AM UTC+1, c.beta...@fz-juelich.de wrote:
>> Dear Helder,
>>
>> OK. Still, there remain two problems: I can choose the color for water and land areas. I did not find any option or keyword to access a certain country. The second problem is, this is not an array that I could modify or fold with anything.
>>
>> Thank you,
>> Clara
>
> Hi,
> has I said, countries or imagery of the earth surface and stuff is not my stuff.
> But I have some good news: here you can find info about how to get coordinates of territories out of shapefile:
> http://www.idlcoyote.com/code_tips/extractpoly.php
> By the way, when you draw the maps of Europe, IDL gets the shapefile (where the borders are defined) from this file:
> IDL path\resource\maps\shape\cntry08.shp
>
> Now the bad news. The pro cgExtractShape is not for free:
> http://www.idlcoyote.com/idldoc/forsale/index.html
>
> You should ask David for a price, might be not that bad... But I sincerely have no idea.
>
> I bought some time ago activecontour.pro and I'm very happy with it.
>
> Cheers,
> Helder
Hi,
if you really *dare* to do the dirty work yourself, here is something to get you started.
Look at the shapefile format here: https://www.esri.com/library/whitepapers/pdfs/shapefile.pdf
Now you can start unveiling what's inside. Here is the header. The rest is left as an exercise for the reader :-)
fn = 'fix_your_path\resource\maps\shape\cntry08.shp'
get_lun, fileUnit
openr, fileUnit, fn
fh = bytarr(100)
readu, fileUnit, fh
free_lun, fileUnit
print, 'File code = ', strtrim(swap_endian(fix(fh[0:3],0,type=13)),2)
print, 'File length = ', strtrim(swap_endian(fix(fh[24:27],0,type=13)),2)
print, 'Version = ', strtrim(fix(fh[28:31],0,type=13),2)
print, 'Shape type = ', strtrim(fix(fh[32:35],0,type=13),2)
print, 'File code = ', strtrim(swap_endian(fix(fh[0:3],0,type=13)),2)
print, 'File length = ', strtrim(swap_endian(fix(fh[24:27],0,type=13)),2)
print, 'Version = ', strtrim(fix(fh[28:31],0,type=13),2)
print, 'Shape type = ', strtrim(fix(fh[32:35],0,type=13),2)
print, 'bounding box xMin = ', strtrim(fix(fh[36:43],0,type=5),2)
print, 'bounding box yMin = ', strtrim(fix(fh[44:51],0,type=5),2)
print, 'bounding box xMax = ', strtrim(fix(fh[52:59],0,type=5),2)
print, 'bounding box yMax = ', strtrim(fix(fh[60:67],0,type=5),2)
print, 'bounding box zMin = ', strtrim(fix(fh[68:75],0,type=5),2)
print, 'bounding box zMax = ', strtrim(fix(fh[76:83],0,type=5),2)
print, 'bounding box mMin = ', strtrim(fix(fh[84:91],0,type=5),2)
print, 'bounding box mMax = ', strtrim(fix(fh[92:99],0,type=5),2)
This is what I got:
File code = 9994
File length = 1906850
Version = 1000
Shape type = 5
bounding box xMin = -180.00000
bounding box yMin = -90.000000
bounding box xMax = 180.00000
bounding box yMax = 83.623600
bounding box zMin = 0.00000000
bounding box zMax = 0.00000000
bounding box mMin = 0.00000000
bounding box mMax = 0.00000000
cheers,
Helder
|
|
|
|
|