map glitch

[mathias]

Hello,

I see a glitch in a map I imported (proc mapimport) from a shapefile.

(see capture.png)

Any idea what causes this and how to correct it ?

The order of the map nodes maybe ?

Here's how I created this map :

map_statisticalsectors weighs 81Mb so I will make a proc gremove on the whole country before.

data map_statisticalsectors;

     set map_statisticalsectors;

     country=1;

     run;

proc gremove data=map_statisticalsectors

     out=map_country all;

     by country;

     id _ID_;

     run;quit;

proc gmap data=map_country

     map=map_country

     all;

     id country;

     choro country / nolegend levels=1 ;

     run;quit;

(see capture2.png)

Thank you,

[GraphGuy]

Hmm ... hard to say...

What's the real id in this map, and are you using an 'id' statement on your Proc Mapimport? The proc will do the best it can without that, but by specifying the id variable you remove a lot of uncertainty & guesswork, which makes a better map. In the example you're setting country=1, but that's not always a valid thing to do, because then gmap might try to inter-connect pieces of the map in the wrong order (thinking they're all just segments of the same id).

For example, here's what happens if you try to plot maps.uscounty similarly, with an invalid id (just making country=1 for all the ids, similar to what you've done in your map)...

data my_map; set maps.uscounty;
country=1;
run;

proc gmap data=my_map map=my_map;
id country;
choro country;
run;

[mathias]

> What's the real id in this map, and are you using an 'id' statement on your Proc Mapimport? The proc will do the best it can without that, but by specifying the id variable you remove a lot of uncertainty & guesswork, which makes a better map. In the example you're setting country=1, but that's not always a valid thing to do, because then gmap might try to inter-connect pieces of the map in the wrong order (thinking they're all just segments of the same id).

The real ID of my map was OBJECTID which I renamed _ID_ .

OBJECTID has the same number of values as adding _ID_ created with CREATE_ID_, so I thought it was ok

Here is my mapimport:

proc mapimport datafile="\\sas9xbi\Data2\ID\DATA\Geo_Spatial\statisticalsectors_belgium\SCBEL01Z4.shp"

  out=imported /*CREATE_ID_*/ ;

  select X Y Area Area_ha CS102001 Nuts3 OBJECTID SUM_PERIME;

  rename OBJECTID=_ID_ CS102001=NIS9 SUM_PERIME=Perimeter;

  run;quit;

adding   "id OBJECTID;" does not change anything visible.

> For example, here's what happens if you try to plot maps.uscounty similarly, with an invalid id (just making country=1 for all the ids, similar to what you've done in your map)...

In my example I do a proc gremove before proc gmap.

Like here : 25591 - Create an outline map of Africa using PROC GREMOVE

data my_map; set maps.uscounty;

    country=1;

    run;

proc gremove data=my_map out=my_map;

    by country;

    id county;

    run;quit;

proc gmap data=my_map map=my_map;

    id country;

    choro country;

    run;

> I narrowed one of the statistical sectors where the first bug appears in my map :

Looks strange...

[GraphGuy]

I think this is still a problem of improper id. In the case of maps.uscounty (which is sort of a special case, because each county requires a compound id of state & county to identify it), the following code would be required to work correctly:

data my_map; set maps.uscounty;
country=1;
run;

proc gremove data=my_map out=my_map;
by country;
id country state county;
run;quit;

proc gmap data=my_map map=my_map;
id country;
choro country;
run;

[GraphGuy]

So, assuming 'objectid' is a good/valid id in your map, I think the following should do what you're wanting ...

proc mapimport
datafile="\\sas9xbi\Data2\ID\DATA\Geo_Spatial\statisticalsectors_belgium\SCBEL01Z4.shp"
out=imported;
id OBJECTID;
run;

data my_map; set imported;
country=1;
run;

proc gremove data=my_map out=my_map;
by country;
id county objectid;
run;

proc gmap data=my_map map=my_map;
id country;
choro country;
run;

[mathias]

This produced the same map as capture2.png in my first post.

I also tried to create a new id with the import :

proc mapimport

datafile="\\sas9xbi\Data2\ID\DATA\Geo_Spatial\statisticalsectors_belgium\SCBEL01Z4.shp"

out=imported create_id_ contents;

select Area CS102001 Nuts3 OBJECTID SUM_PERIME;

*id OBJECTID;

run;

data my_map; set imported;

country=1;

run;

proc gremove data=my_map out=my_map;

by country;

id country _ID_;

run;

proc gmap data=my_map map=my_map;

id country;

choro country;

run;

and to add more ID's :

id OBJECTID CS102001 Nuts3 Sector_fr Sectie_fr SEC102001;

but it always result in capture2.png

[GraphGuy]

If there's some way you could get the shapefiles to me, I'd be happy to experiment with them and see if I can read them in. Is there a website (or ftp site) I can download them from? (sounds like they're too big to email). If you don't want to divulge the location to the whole group, feel free to email it to me at robert.allison@sas.com

[GraphGuy · Accepted solution]

Ok - I've done a little experimenting with your shapefile, and I believe those aren't "glitches" but actually how the map is supposed to look

For example, that odd little rope-like area that goes off to the east side of Belgium, when I look up that area in Google maps, it shows that to be a road, that Belgium apparently owns the land to the immediate left & right of. Zoom into the following to see:

belgium country - Google Maps

I had a friend use some non-SAS software to view the ESRI shapefiles directly, and here is a screen-capture of what they look like:

I imported it into SAS, using Proc Mapimport, but with 'id nuts1 nuts2 nuts3' (rather than having mapimport create an id for each area), and the gmap of this one looks like your gmap (and like the non-SAS map).

I then used Proc Gremove to remove all the internal boundaries, and got the following (which looks similar to yours):

Here's the code I used, in case you want to use the same:

%let name=map;
filename odsout '.';

libname foo 'c:\public\tmp\';

proc mapimport datafile="./SCBEL01Z4.shp" out=foo.my_map;
id Nuts1 Nuts2 Nuts3;
run;

data my_map; set foo.my_map (keep = X Y Nuts1 Nuts2 Nuts3 SEGMENT);
run;

goptions device=png;

ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" (title="Map") style=htmlblue;

title "Original map, with all variables";
proc print data=foo.my_map (obs=10);
run;

title "Original Map, with 'id nuts1 nuts2 nuts3'";
proc gmap data=my_map map=my_map;
id nuts1 nuts2 nuts3;
choro segment / levels=1 nolegend
des='' name="&name";
run;

data removed; set my_map;
country=1;
run;

proc gremove data=removed out=removed;
by country;
id country nuts1 nuts2 nuts3;
run;

title "Added country=1, and then gremoved internal boundaries";
proc gmap data=removed map=removed;
id country;
choro country / levels=1 nolegend
des='' name="&name";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;

[GraphGuy]

The last section in this chapter has some tips on using Proc Mapimport that might be helpful:

http://robslink.com/SAS/book1/Chapter_05_Maps.pdf