cogo - Data format for checking a land survey deed

A commercial property is in the closing phase of sale and the only hold up is that the deed description is one that's been lifted from prior deeds dating back to the 1950's and is obviously erroneous. A 1992 survey does appear to properly depict what the Seller's have to convey and both Buyers and Seller agree that the aforementioned survey be the basis of description in the new deed and title policy if ...it plots out in reasonable fashion in Copan.

But, after spending six hours trying to get a result (closure error, or an additional call necessary to close) and a graph, I have been unsuccessful. format I've used is:

Map Traverse
1 degrees distance tab tab tab
2 degrees distance tab tab tab
3 degrees distance tab tab tab
1 degrees distance tab tab tab

That produced a bizarre graph, so I tried making the first call

Map Traverse

1 0.0 0.0 tab tab tab

Then I get an error message saying in effect "insufficent tabs" and use cut and paste to provide missing tabs.

---

My data:

1                   
2   354.0847    129.65          
3    84.4149    117.82          
4   352.2410     69.78          
5   309.4352     82         

6    37.1300     55         
7   307.4700     61         
1   276.5342     95.62          

---

The original survey plat:

arcgis desktop - Find out which raster has highest values

I have 3 rasters for water use and 1 raster for water supply.

I want to know which one has the highest values which will mean that it has the greatest influence on available water supply.

What tools can I use for this in arcmap 10.0 and how do I do this?

ESRI REST API clients?

So, if I implement the esri GeoServices REST API with a non-esri service, what clients can I use with it?

We (GeoREST) were initially thinking that ArcGIS Explorer or other desktop clients would be able to access it, but this doesn't seem to be the case. Looking at the web stream with Fiddler, these appear to fall back to a SOAP interface when passed the REST end point.

The web APIs looked promising, but the terms under which you're allowed to use them appear to be pretty restrictive:

I thought these permitted use for public, non-commercial purposes, but wording in the second link made me doubt that interpretation.

I may be missing something, but a lack of freely usable clients seems to limit the value for a third party to implement this API.

Does anyone have any further insight into the licencing restrictions on the web APIs, or know of other services that grok the REST API?

carto - Getting info from multiple rows in single popup of CartoDB?

If you have a table like this:

Name | Item | State ---------|--------|----------- Alice | Apple | Alabama Ben | Banana | Alabama Caitlin | Cherry | California

Assuming you have created polygons based on the states, how could you create an 'info window' which combines the details of all rows in the same state?

So for example, the popup for Alabama might be something like:

State: Alabama
People: 2
Items: Apple, Banana

As far as I can tell, the default behaviour for Cartodb would be to place one polygon on top of the other, and for the info window to display the info from the top most polygon only.

openlayers - Save features and attributes as WFS in GeoServer

I want to create some squared features with OpenLayers 4, add attributes to each feature and want to save them on the GeoServer to load them later on to another map. I am able to create the features, add attributes to them, read the attributes again and load a layers as WFS again.

My only problem is to save the created features to the GeoServer. I figured out that saving and loading should be possible with WFS-t. I found some examples online but they don't seem to work properly.

When i send the feature to the GeoServer with the function transactWFS at line 133 I get the Error "org.geoserver.wfs.WFSTransactionException: Feature type 'Stairways' is not available" in the GeoServer log file. The Layer Stairways exists on the GeoServer. This is the XML file that is sent:

  
    
      
        
          
            
              1528915.8154546698
                    6627851.731194374 1528917.9110709063

                    6627858.932471524 1528925.1123480562
                    6627856.836855288 1528923.0167318196
                    6627849.635578138 1528915.8154546698
                    6627851.731194374
            
                                                             
        
      
    
         

This is my whole Code:

   GridSource.addFeature(NewSquare);

    var attr = {};
    var key = "TestAttr";
    attr[key] = "This is a Test";

    NewSquare.attributes = attr; 

    transactWFS('insert', NewSquare);        
  }



  var formatWFS = new ol.format.WFS();

  var formatGML = new ol.format.GML({
      featureNS: gs.wfs,
      featurePrefix: 'HTW_Erd',

      featureType: 'Stairways',
      srsName: 'EPSG:3857'
  });

  var xs = new XMLSerializer();


  var transactWFS = function (mode, f) {
      var node;
      switch (mode) {

          case 'insert':
              node = formatWFS.writeTransaction([f], null, null, formatGML);
              break;
          case 'update':
              node = formatWFS.writeTransaction(null, [f], null, formatGML);
              break;
          case 'delete':
              node = formatWFS.writeTransaction(null, null, [f], formatGML);
              break;
      }

      var payload = xs.serializeToString(node);
      $.ajax(gs.ows, {
          service: 'WFS',
          type: 'POST',
          dataType: 'xml',
          processData: false,
          contentType: 'text/xml',
          data: payload,
          error: function(e) {
              var errorMsg = e? (e.status + ' ' + e.statusText) : "";

              alert('Error saving this feature to GeoServer. \n\n'
                  + errorMsg);
          }
      }).done();

  };

Edit:
I saw in the log file that i wasnt allowed to write on that layer so i set the permissions. I also changed the projection of the GridLayer to 3857 and the name of the feature to geometry. But that didn't help either. I now don't get any errors and the log looks like this:

2017-07-30 23:47:07,676 INFO [geoserver.wfs] - 

Request: getServiceInfo
2017-07-30 23:47:07,709 INFO [geoserver.wfs] - 
Request: getServiceInfo
2017-07-30 23:47:07,744 INFO [geoserver.gwc] - DataStoreChange: {https://localhost:8080/geoserver/htw_erd}Stairways PreInsert
2017-07-30 23:47:07,749 INFO [geoserver.gwc] - DataStoreChange: {https://localhost:8080/geoserver/htw_erd}Stairways PostInsert
2017-07-30 23:47:07,751 INFO [geoserver.wfs] - 
Request: getFeature
    service = WFS
    version = 1.1.0
    baseUrl = http://localhost:8080/geoserver/

    query[0]:
        filter = [  bbox POLYGON ((1528628.1984057308 6627469.038650171, 1528628.1984057308 6628141.798427672, 1529548.9368749668 6628141.798427672, 1529548.9368749668 6627469.038650171, 1528628.1984057308 6627469.038650171)) ]
        srsName = EPSG:3857
        typeName[0] = {https://localhost:8080/geoserver/htw_erd}Stairways
    outputFormat = application/json
    resultType = results
2017-07-30 23:47:07,754 INFO [wfs.json] - about to encode JSON
2017-07-30 23:47:07,801 INFO [geoserver.wfs] - 
Request: transaction
    service = WFS

    version = 1.1.0
    baseUrl = http://localhost:8080/geoserver/
    group[0] = wfs:insert=net.opengis.wfs.impl.InsertElementTypeImpl@e314ac (feature: [SimpleFeatureImpl:Stairways=[SimpleFeatureImpl.Attribute: the_geom=null, SimpleFeatureImpl.Attribute: id=null, SimpleFeatureImpl.Attribute: indoor_are=null, SimpleFeatureImpl.Attribute: indoor_lev=null, SimpleFeatureImpl.Attribute: indoor_l_1=null]], handle: null, idgen: , inputFormat: , srsName: null)
    insert[0]:
        feature[0] = SimpleFeatureImpl:Stairways=[SimpleFeatureImpl.Attribute: the_geom=null, SimpleFeatureImpl.Attribute: id=null, SimpleFeatureImpl.Attribute: indoor_are=null, SimpleFeatureImpl.Attribute: indoor_lev=null, SimpleFeatureImpl.Attribute: indoor_l_1=null]
        idgen = GenerateNew
        inputFormat = text/xml; subtype=gml/3.1.1
    releaseAction = ALL

All the features that are already in the Layer named 'the_geom', renaming the feature to the_geom doesn't help.

Edit 2: This is the Output

No service: ( wfs\ )

I have to add that I dont have the default layers anymore, so i dont have the topp workgroup and the layer states anymore.

Edit 3: This is the Output for this link: http://localhost:8080/geoserver/wfs?service=wfs&version=1.1.0&request=DescribeFeatureType&typeName=HTW_Erd:Stairways

Answer

You are telling GeoServer that the namespace for your feature is xmlns="http://localhost:8080/geoserver/wfs" (or gs.wfs) when it should be the URI associated with your workspace that contains the Stairways featuretype.

I would also check to be sure that the layer name really starts with a capital S as it is case sensitive.

EDIT

You will need to set the geometry field name using

 geometryName: "the_geom",

in your GML format.

arcpy - Get a point on a polyline from end points given distance along poly line

I have start point and end point of a polyline. How can I get a point on that polyline from end points specified by the given distance.

How can I get this by using arcpy provided that this script should work on ArcView without importing any modules (Like this answer). I'm trying not to hard code values.

Related Question:

Seeking point cloud (LiDAR) data?

Is there any freely available pointcloud data?

I'm especially interested in dense clouds of small areas(5m to 100m radius) but anything else would be fine too.

I've found that these pages have very nice and easy to access point cloud data.

Answer

I decided to merge other answers with mine and organize them into a tabular format. I think it is easier to read and manage for future visitors:

The table can be accessed from the following link in csv format:

View in tabular form: Free LiDar DataSources

Download (csv): Free LiDar DataSources

Please submit pull request if you intended to add to this list.

sld - Mixed geometry types styling with GeoServer 2.2

I have a PostGIS table with geometry column. In this column I have geometries of several types (polygons, points, lines). I publish this this table on GeoServer.

Now I want to create OpenLayers layer with this data. But how to make a style for this data?

I've already read http://docs.geoserver.org/stable/en/user/styling/sld-tipstricks/mixed-geometries.html and don't want to create gtype in PostGIS if this possible.

I try to create new style within GeoServer but I think that do something wrong:

 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"

xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

    
     
      
       the_geom
      
      Polygon
     

     
       
          #0000ff
          1
        
     
    

And get a error:

line 10: cvc-complex-type.2.4.a: Invalid content was found starting with element 'Rule'. One of '{"http://www.opengis.net/sld":Name, "http://www.opengis.net/sld":Title, "http://www.opengis.net/sld":Abstract, "http://www.opengis.net/sld":NamedLayer, "http://www.opengis.net/sld":UserLayer}' is expected.
line 11: cvc-complex-type.2.4.a: Invalid content was found starting with element 'ogc:PropertyIsEqualTo'. One of '{"http://www.opengis.net/sld":Name, "http://www.opengis.net/sld":Title, "http://www.opengis.net/sld":Abstract, "http://www.opengis.net/sld":LegendGraphic, "http://www.opengis.net/ogc":Filter, "http://www.opengis.net/sld":ElseFilter, "http://www.opengis.net/sld":MinScaleDenominator, "http://www.opengis.net/sld":MaxScaleDenominator, "http://www.opengis.net/sld":Symbolizer}' is expected.

What am I doing wrong?

UPDATE

I try to show a points and make this rule for this geom type:

 
                
            
                

                    the_geom
                
                Point
            
            
        
        
            #0000ff
            2
        

        
    

But not see then on my map. What's wrong?

Answer

The first validation error is because you are missing Name, Title, Abstract, NamedLayer and UserLayer elements in your SLD. The second validation error refers to the missing filter element inside your Rule element. So if you add the filter element and the standard elements in the sld elemet your SLD should validate fine. Finally our SLD should look like this:

  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">

  
    
    
      
      
      
        
                    
                
                    

                        the_geom
                    
                    Polygon
                
                
            
            
                #0000ff
                1
            

            
        
      
    
  

If you want to style other geometrytypes - just create a new rule with the symbolizer and geometrytype filter you want to style.

UPDATE

Your point rule should have a PointSymbolizer element instead of the PolygonSymbolizer. So your point rule look like this:

 
                
            
                
                    the_geom
                
                Point
            
            
        

        
              
                circle
                
                  #0000ff
                  2
                
              
              6
            

        
    

See more on how to style points in the GeoServer SLD Cookbook

dem - Viewing LIDAR data (.las) in QGIS?

Is there an easy path to visualizing LIDAR data in QGIS?

I have some USGS LIDAR data in .las format downloaded from http://lidar.cr.usgs.gov/. This means I have both the .las and metadata in .xml format. I am aware liblas, but not how to apply it to this task. I am running on Ubuntu 11.04 with QGIS 1.7.0-Wroclaw.

A similar question for ArcGIS is: Converting LiDAR data to raster (DEM/DSM) for ArcGIS input?. I just need to get a sense of this data and the registration. Conversion to a DEM would be OK if I can visualize it.

Performing multiple layer intersection in QGIS?

I am trying to find a way with QGIS to get the intersection area of more than 2 layers.

All the tools in QGIS take just the 2 input layers.

For example in the image below, I have 4 layers and I want to retrieve just the black area.

Furthermore, I need to retrieve the count attribute and add it to the intersected area attribute table (in this case 4)

I know there is a tool that does the trick in ArcGIS, but is there a QGIS solution?

Answer

Forgot to post my solution...
STEP 1 : Merge layers
In the end I changed my Qgis Version to 2.14.11, the algorithm takes a single input param (layers) separated by a semicolon :
processing.alghelp('qgis:mergevectorlayers')

ALGORITHM: Merge vector layers
    LAYERS 
    OUTPUT 

Prior it was :

processing.runalg('qgis:mergevectorlayers', layer1, layer2, output)

You can also do a loop to use this version of the algorithm

for i, layer in enumerate(layers):

            if i == 0:
                layer1 = layers[i]
                layer2 = slayers[i+1]
                processing.runalg('qgis:mergevectorlayers', layer1, layer2, savename)
            if i == leng-1:
                break

            else:
                layer1 = QgsVectorLayer(savename, 'temp', 'ogr')
                layer2 = sme_layers[i+1]
                savename = name          

                processing.runalg('qgis:mergevectorlayers', layer1, layer2,  savename)

STEP2 : overlay the layer with itself

processing.alghelp('grass:v.overlay')




ALGORITHM: v.overlay - Overlays two vector maps.
    ainput 
    atype 
    binput 
    operator 
    -t 
    GRASS_REGION_PARAMETER 
    GRASS_SNAP_TOLERANCE_PARAMETER 

    GRASS_MIN_AREA_PARAMETER 
    GRASS_OUTPUT_TYPE_PARAMETER 
    output 

STEP 3 Joseph answer seems decent

batch - How to use QGIS Atlas with multiple layers/files per coverage?

I've got 30 contour datasets that all sit over the same extent. Each contour set is in its own file. I also have 30 extent polygons for each contour set as a separate file.

Can I use Atlas to generate a series of maps - each displaying a 'pair' of datasets - ie. contours and extent for each pair over the same extent?

Everything I've read indicates that you can loop through features in a layer with Atlas, but not through multiple layers. Could I work around this by adding a 'layer display list' as an attribute of the coverage layer?

Maybe with a python script?

---

Update:

This is going back to a previous question from a few years ago that I asked and didn't find a solution for. QGIS has advanced a few versions since then - so I'm hoping there is either an updated solution, or someone could provide guidance to develop a script.

Advanced warning - I'm a geologist - layers of rock are how I think! :-)

I have a project area for which I want to produce a series of maps (an atlas). Each map will cover exactly the same extent, but will show properties for a different layer of rock. In this example - contours of the depth to the top of 50 different rock layers.

The desired output would be 50 maps, same extent, different layer (or layer set) displayed

To add more complexity - I'd like to add 'sets' of layers. Rock layer 1 also has unique shapefiles showing the extent polygon layer and also a label_point layer.

Current Method

I've tried merging all the contour layers into one shapefile and using this as the Atlas coverage layer. This method doesn't quite produce the desired result because:

• I can't filter the combined layer to only show one layer (not one feature) at a time.



• I can filter using this rule based style expression on the combined layer! (where unit is the rock layer name):

attribute($currentfeature , 'unit') = attribute(@atlas_feature, 'unit')

• Critical issue - I get one page per individual contour (ie. one per feature - so hundreds of pages per rock layer). I can't work out how to filter so I only get one page per rock layer (or if that is possible).


• I can't merge the different 'sets' (points, polylines, polygons). Workaround - I could merge and filter based on geometry type - same as above.


• Merging 50 input files into one master shapefile seems a bit redundant - surely one more code loop could negate the need for this? Workaround - use the processing framework to combine input files.

Answer

Since QGIS 2.12, with the introduction of data-defined properties on the "lock layers for map" option, in the map item properties, you can do what you seek. You just need to create the right coverage layer, with the information of the extents, but also of the layers to "print".

Create the layer presets

You can "save" the current state (visible or not visible) of your layers as presets. Just create the desired combination and press the eye icon on the layers panel, and choose "add preset".

Basically, you will have to create a preset for each contour, with that layer visible (paired with other layers if you like). This will be your "set of layers".

Give your preset memorable names, or else just opt for some kind of a sequence.

Creating the coverage layer.

Now you will have to create a polygon layer with all the possible combinations between the extents polygons, and all the presets. That means you will end up with some redundancy. You will have to replicate each extent polygon for each preset. Something like this:

polygon 1, preset a polygon 1, preset b ... polygon 50, preset a polygon 50, preset b

You can do this manually with copy paste, or else create some kind of script, but I would say your best and easier way would be using Spatialite or Postgis database. Just import you polygons extents there, and create a new non-spatial table with all presets names. Then, create a view that returns all the combinations.

SELECT f.*, g.preset_name
FROM polygons_extent as f, presets as g

Making it work

Add the new layer or view to QGIS, use it as coverage layer in atlas settings, and use the preset_name attribute in data-defined properties on the "lock layers for map" option.

qgis - What's the purpose of .qgs~ files?

When i work on QGIS project (named "pro"), some file with this symbol appears in the workspace folder with the name "pro.qgs~":

I tried to open it- but with no success.

1. How can i open it and what is the purpose of it?


2. What to do, so this file will not create in the folder?

Answer

pro.qgs~ is a backup copy of pro.qgs. See: When saving a .qgs file in QGIS, another appears. Is this a backup?

If your project file gets corrupted for any reason, you can simply rename pro.qgs~ to pro.qgs and continue to work from there.

Like the original .qgs file, you can open it with a text editor.

Why would you want to prevent creation of a backup file?

Measuring accuracy of latitude and longitude?

I have latitude and longitude as 19.0649070739746 and 73.1308670043945 respectively.

In this case both coordinates are 13 decimal places long, but sometimes I also get coordinates which are 6 decimal places long.

Do fewer decimal points affect accuracy, and what does every digit after the decimal place signify?

Answer

Accuracy is the tendency of your measurements to agree with the true values. Precision is the degree to which your measurements pin down an actual value. The question is about an interplay of accuracy and precision.

As a general principle, you don't need much more precision in recording your measurements than there is accuracy built into them. Using too much precision can mislead people into believing the accuracy is greater than it really is.

Generally, when you degrade precision--that is, use fewer decimal places--you can lose some accuracy. But how much? It's good to know that the meter was originally defined (by the French, around the time of their revolution when they were throwing out the old systems and zealously replacing them by new ones) so that ten million of them would take you from the equator to a pole. That's 90 degrees, so one degree of latitude covers about 10^7/90 = 111,111 meters. ("About," because the meter's length has changed a little bit in the meantime. But that doesn't matter.) Furthermore, a degree of longitude (east-west) is about the same or less in length than a degree of latitude, because the circles of latitude shrink down to the earth's axis as we move from the equator towards either pole. Therefore, it's always safe to figure that the sixth decimal place in one decimal degree has 111,111/10^6 = about 1/9 meter = about 4 inches of precision.

Accordingly, if your accuracy needs are, say, give or take 10 meters, than 1/9 meter is nothing: you lose essentially no accuracy by using six decimal places. If your accuracy need is sub-centimeter, then you need at least seven and probably eight decimal places, but more will do you little good.

Thirteen decimal places will pin down the location to 111,111/10^13 = about 1 angstrom, around half the thickness of a small atom.

Using these ideas we can construct a table of what each digit in a decimal degree signifies:

• The sign tells us whether we are north or south, east or west on the globe.


• A nonzero hundreds digit tells us we're using longitude, not latitude!


• The tens digit gives a position to about 1,000 kilometers. It gives us useful information about what continent or ocean we are on.


• The units digit (one decimal degree) gives a position up to 111 kilometers (60 nautical miles, about 69 miles). It can tell us roughly what large state or country we are in.


• The first decimal place is worth up to 11.1 km: it can distinguish the position of one large city from a neighboring large city.


• The second decimal place is worth up to 1.1 km: it can separate one village from the next.


• The third decimal place is worth up to 110 m: it can identify a large agricultural field or institutional campus.


• The fourth decimal place is worth up to 11 m: it can identify a parcel of land. It is comparable to the typical accuracy of an uncorrected GPS unit with no interference.


• The fifth decimal place is worth up to 1.1 m: it distinguish trees from each other. Accuracy to this level with commercial GPS units can only be achieved with differential correction.


• The sixth decimal place is worth up to 0.11 m: you can use this for laying out structures in detail, for designing landscapes, building roads. It should be more than good enough for tracking movements of glaciers and rivers. This can be achieved by taking painstaking measures with GPS, such as differentially corrected GPS.



• The seventh decimal place is worth up to 11 mm: this is good for much surveying and is near the limit of what GPS-based techniques can achieve.


• The eighth decimal place is worth up to 1.1 mm: this is good for charting motions of tectonic plates and movements of volcanoes. Permanent, corrected, constantly-running GPS base stations might be able to achieve this level of accuracy.


• The ninth decimal place is worth up to 110 microns: we are getting into the range of microscopy. For almost any conceivable application with earth positions, this is overkill and will be more precise than the accuracy of any surveying device.


• Ten or more decimal places indicates a computer or calculator was used and that no attention was paid to the fact that the extra decimals are useless. Be careful, because unless you are the one reading these numbers off the device, this can indicate low quality processing!

postgis - pgr_dijkstra gives wacky routes sometimes with undirected graph

I am using the function below to find the shortest path between several origin-destination pairs, and it seems to work correctly most of the time.

CREATE OR REPLACE FUNCTION public.sp_od(
    orig integer,
    dest integer)
    RETURNS TABLE(shortest_path geometry) 
LANGUAGE 'sql'


AS $BODY$   

SELECT st_makeline(geom) as shortest_path
FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom, true) as cost FROM public."WA_roads"',
(SELECT source FROM public."WA_roads"
 ORDER BY ST_StartPoint(geom) <->
     (select ST_SetSRID(ST_MakePoint(CAST(ocentx as double precision), CAST(ocenty as double precision)), 4326) from all_trips_non_zero where origin = orig LIMIT 1) ASC
LIMIT 1),

(SELECT source FROM public."WA_roads"
 ORDER BY ST_StartPoint(geom) <->
     (select ST_SetSRID(ST_MakePoint(CAST(dcentx as double precision), CAST(dcenty as double precision)), 4326) from all_trips_non_zero where destination = dest LIMIT 1) ASC
LIMIT 1), directed := false
) as pt
JOIN public."WA_roads" rd ON pt.edge = rd.id;

$BODY$;

The above function takes an origin and destination zip code and using the corresponding lat, long for the zip code from 'all_trips_non_zero' table, finds the nearest "source" node in our road network closest to the origin and destination points to use for shortest path calculation.

One of the problematic path is shown in the image below. Origin in red and destination is yellow. Brown is the network and blue is the shortest path from the above function.

For reference, the shortest path predicted by Google is as under.

I cannot explain some of the straight lines, not on the network, in the pgr_dijkstra path above. These are somehow part of the geometry and I get a total path length of around 360 miles, which is more than twice the distance predicted by Google. When I use directed:= true, I get no path for this OD pair. The "WA_roads" shapefile and "all_trips_non_zero" file are here. I further transformed the shapefile to EPSG:4326 as described here. The pgr_analyzegraph output is below:

NOTICE:  pgr_analyzeGraph('WA_roads',1e-06,'geom','id','source','target','true')
NOTICE:  Performing checks, please wait ...
NOTICE:  Analyzing for dead ends. Please wait...
NOTICE:  Analyzing for gaps. Please wait...
NOTICE:  Analyzing for isolated edges. Please wait...

NOTICE:  Analyzing for ring geometries. Please wait...
NOTICE:  Analyzing for intersections. Please wait...
NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 0
NOTICE:                            Dead ends: 214
NOTICE:  Potential gaps found near dead ends: 0
NOTICE:               Intersections detected: 1
NOTICE:                      Ring geometries: 2

Successfully run. Total query runtime: 1 secs 204 msec.

1 rows affected.

Raster calculator reclass 0 as 0 QGIS

I want to process a raster with values 0-16. I want to set all values between 1 and 15 as 1, a value of 16 as 0 and cells with a value of 0 as 0.

How do I maintain a cell value of 0?

I have tried the following:

("raster_name@1" !=0 AND "raster_name@1" <= 15*1)+("raster_name@1" = 16*0)

Results in: Values of 16 as 0, values 1-15 as 1, values 0 also = 1

I also tried:

("raster_name@1" !=0 AND "raster_name@1" <= 15*1)+("raster_name@1" 16*0)+("raster_name@1" = 0*0)

Results in: Values of 16 as 0, values 1-15 as 1, values 0 = 2

Please, can you explain how to retain cells with a 0 value as 0?

Answer

Simply write "raster_name@1" != 0 AND "raster_name@1" <= 15 in order to set all values between 1 and 15 to 1, and everything else to 0.

Exporting specified layers in ArcMap using ArcPy?

I am trying to export specific layers from an ArcMap .mxd. From the TOC below I would like to export "irrig_2, citylyr, boundlyr" as a GIF. Once this is done the script would then export "irrig_1, citylyr, boundlyr" as another GIF and so on. If possible, the script would also export a legend that corresponds to the layers being exported. I'm not sure if it's a good place to start, but I'm trying to work off the code below.

import arcpy

mxd = arcpy.mapping.MapDocument(r"E:\IrrigatedLands\FC_test\irrig2.mxd")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]

GIFPath = r"E:\IrrigatedLands\FC_test"

for lyr in arcpy.mapping.ListLayers(mxd, '', df):
        if lyr.name == "citylyr":

            lyr.visible = True
        if lyr.name == "boundlyr":
            lyr.visible = True
        if lyr.name == lyr:
            lyr.visible = True
            arcpy.mapping.ExportToGIF(mxd,GIFPath+"\\" + lyr.name + ".gif")
            lyr.visible = False
arcpy.RefreshActiveView()
del mxd

Answer

You need to loop through a list of your output layers so you can turn them on and off for each export.

import arcpy

projFolder = r"E:\IrrigatedLands\FC_test"

mxd = arcpy.mapping.MapDocument(r"{}\irrig2.mxd".format(projFolder))
df = arcpy.mapping.ListDataFrames(mxd, '')[0]

GIFPath = r"E:\IrrigatedLands\FC_test"


layers = ['irrig_0','irrig_1','irrig_2'] # List of your output layers

lyrlist = arcpy.mapping.ListLayers(mxd, '', df)

for layer in layers:
    layerOn = "" # When required layer is turned on this will be the layer name to output
    for lyr in lyrlist:
        lyr.visible = False
        #if lyr.name in layers and lyr.name <> layer:
        #   lyr.visible = False

        if lyr.name == "citylyr":
            lyr.visible = True
        if lyr.name == "boundlyr":
            lyr.visible = True
        if lyr.name == layer:
            lyr.visible = True
            layerOn = lyr.name
    if layerOn:
        arcpy.mapping.ExportToGIF(mxd,r"{}\{}.gif".format(GIFPath,layerOn))


del mxd

carto - CartoDB only georeferences half of my table

I imported a CSV file into CartoDB, with over a million points, all with a clear (and correct) latitude and longitude. For many of the collected data points the georeference column says ‘null’, while the latitude & longitude are correct.

So CartoDB only georeferenced half of my table and I am missing a lot of points.

I do have a payed account and did not reach my maximum of points I can georeference.

Anybody know how I can get my whole CSV georeferenced?

arcgis desktop - How to get values of each cell in raster attribute table?

How do I get the values of each cell in the raster attribute table. Usually the attribute table gives you the value and number of counts. But, in my case I want each cell listed and corresponding cell values.

Random sampling of raster using R?

Is there a straightforward way of randomly sampling a raster so that the output of the process is a raster?

I'm using an example that I found on the r-sig-geo list and I have also tried the sampleRandom function in the raster package. Both of these approaches produce an output that I am not certain how to transform into a raster. I was not able to find an approach after searching for several combinations of "SpatialPointsDataFrame raster".

library(raster)

# read in raster
rasterSource <- 'landsat.TIF'

r <- raster(rasterSource)

# convert to spatial points data frame
r.spgrd<-as(r,"SpatialPointsDataFrame") 

# elminate NA values
r.spgrd = r.spgrd[!is.na(r.spgrd[[1]]),] 

# sample points
selectedPoints = sample(1:length(r.spgrd[[1]]), 1000)

r.sampled = r.spgrd[selectedPoints,]

# try to make spgrd into a raster
r.test <- raster(r.sampled)

When I run r.test I get the output:

class       : RasterLayer 
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 28617, 14766  (x, y)
extent      : 1838505, 2124675, 2328685, 2476345  (xmin, xmax, ymin, ymax)

coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
values      : none

So that the following line which tries to write a raster produces the message:

# write out as ascii file
writeRaster(r.test, filename="test1.ASC", datatype="ascii", overwrite=TRUE)

Error: hasValues(x) is not TRUE

My main objective is to produce some type of raster after the sampling process. I'm also fine with just changing the values within my raster (I'm just not certain how to do that).

Answer

You can adapt examples from the Raster package vignette, section 5.2. Here's one way:

r <- raster(ncol=30,nrow=20)
r[] <- 1:(30*20)              # Raster for testing
#plot(r)                      # (If you want to see it)
r[runif(30*20) >= 0.30] <- NA # Randomly *unselect* 70% of the data
plot(r)

python - I have NDVI layer in Earth Engine. I want to assign and calculate area of good and bad ndvi area

    var mergedCollection = s2.map(function(image) {
    return image.select(['B8', 'B4'], ['NIR', 'red'])
    })


   var mergedAllFunction = function(image) {
   var ndvi = image.normalizedDifference(['NIR', 'red']).rename('NDVI');



   var thres1 = ndvi.gte(0).rename('thres1')
   var thres2=ndvi.gt(0.1) && ndvi.lt(0.3).rename('thres2')
   var thres3=ndvi.gt(0.4) && ndvi.lt(0.7).rename('thres3')
   var thres4=ndvi.gt(0.8) && ndvi.lte(1).rename('thres4')

   return image.addBands(ndvi).addBands(thres1,thres2,thres3,thres4);
   print(ndvi)
   }

How do I get area (count of pixels under each threshold)

Answer

There are a couple problems in your code. First, using && in that way will not work for reasons described here. Second, you can't addBands(arg1, arg2,...). Those args should be in a list. Here's a complete example that incorporates your code and a previous answer:

var mergedCollection = s2
    .filterDate('2018-01-01', '2018-03-01')
    .map(function(image) {
      return image.select(['B8', 'B4'], ['NIR', 'red'])
    })

var mergedAllFunction = function(image) {

 var ndvi = image.normalizedDifference(['NIR', 'red']).rename('NDVI');

 var thres1 = ndvi.gte(0).rename('thres1')
 var thres2=ndvi.gt(0.1).and(ndvi.lt(0.3)).rename('thres2')
 var thres3=ndvi.gt(0.4).and(ndvi.lt(0.7)).rename('thres3')
 var thres4=ndvi.gt(0.8).and(ndvi.lte(1)).rename('thres4')

 return image.addBands(ndvi).addBands([thres1,thres2,thres3,thres4]);
}


var median = mergedCollection.median();

var merged = mergedAllFunction(median);

var areas = merged
    .select(['thres1', 'thres2', 'thres3', 'thres4'])
    .multiply(ee.Image.pixelArea())
    .reduceRegion({
      reducer: ee.Reducer.sum(),
      geometry: geometry,  // a geometry

      scale: 10,   // scale = 10 for sentinel-2 'red' band
      maxPixels: 1e9  
    });
print(areas)

arcgis desktop - Using line of sight or something similar to identify where lines go through hills

When QCing power infrastructure I am having to zoom into every gap and check to make sure the line makes sense (not going through the mountain). I typically check for odd spacing between the pylons and zoom into the imagery. What I want to accomplish is have some tool identify where the 'line of sight' would intersect with the DTED_Level_2. Is there a tool in ArcMap 10.1 that can accomplish this? I have 3D Analyst installed.

Enabling CORS in GeoServer (jetty)?

I hope somebody has already figured this one out. I just installed Geoserver 2.9 on a vanilla Ubuntu 16.04 distro. The Geoserver 2.8 method of enabling CORS with the shanbe.hezoun class does no longer work with Jetty 9.2.13.

There are mentions that CORS support is already packaged with Jetty 9.2.13 in the jetty-servlets.jar.

The Jetty lib which is compiled with Geoserver contains a jetty-servlet-9.2.13.v20150730.jar in geoserver/lib but not jetty-servlets.9.2.13.v20150730.jar. Are these supposed to be the same jar with a different name?

It should be possible to enable CORS either in geoserver/etc/webdefault.xml or in geoserver/webapps/geoserver/WEB-INF/web.xml.

My understanding is that the webdefault.xml is applied first and the web.xml thereafter.

I have tried following filter in both xml. I haven't got as far as adding a filter mapping. Adding the filter alone will cause the Geoserver/Jetty service to not start proper.

    cross-origin
    org.eclipse.jetty.servlets.CrossOriginFilter

Answer

Edit the webapps/geoserver/WEB-INF/web.xml file. There are two references to CORS in this file:

  cross-origin
     org.eclipse.jetty.servlets.CrossOriginFilter
  

and

   cross-origin
   /*

You must uncomment both blocks (that is remove from the filter and filter-mapping blocks.

Then when you restart Jetty you can test that everything is working by using a command like:

curl -v -H "Origin: http://example.com" http://astun-desktop:9080/geoserver/wfs\?service\=WFS\&version\=2.0.0\&request\=GetFeature\&typenames\=sf:bugsites\&filter\=%3Cfes:Filter%20xmlns:fes\=%22http://www.opengis.net/fes/2.0%22%3E%3Cfes:ResourceId%20rid\=%22bugsites.3%22/%3E%3C/fes:Filter%3E

which if all is well will give a result like:

> User-Agent: curl/7.35.0

> Host: astun-desktop:9080
> Accept: */*
> Origin: http://example.com
>  
< HTTP/1.1 200 OK 
< Access-Control-Allow-Origin: http://example.com 
< Access-Control-Allow-Credentials: true 
< Access-Control-Expose-Headers:  
< Content-Type: text/xml; subtype=gml/3.2 
< Content-Disposition: inline; filename=geoserver-GetFeature.text 

< Transfer-Encoding: chunked
* Server Jetty(9.2.13.v20150730) is not blacklisted 
< Server: Jetty(9.2.13.v20150730) 
< 
* Connection #0 to host astun-desktop left intact 
590529 49146253Beetle site%

Update 24th Oct 2019

It it is no longer necessary to add the following jar to GeoServer (at least with versions 2.13.x and later) and it will cause an error. I'm leaving this note here for people fighting older versions.

1. Add the Jetty-Utility Servlets Jar to match the version of Jetty - for current versions of GeoServer (2.15.x) it is 9.4.12.v20180830, copy this to webapps/geoserver/WEB-INF/lib inside the geoserver-2.15.0 directory (or wherever you unpacked the zip file).

geotiff tiff - Splitting .tif image into several tiles?

I have an image of size 1GB (.tif), with the width and height 94000x71680. I would like to chunk this image into 20000X20000 tiles so that I can process them.

How can I do this?

Answer

I propose two solutions: the first one using QGIS, the second one using Python (GDAL).

---

Solution using QGIS

In QGIS you may create a VRT mosaic.

Please follow this procedure (see the image below):

1. Load the raster in the Layers Panel;



2. Right-click on it and choose Save As...;


3. Check the Create VRT option;


4. Choose the folder where your outputs will be saved;


5. Set the extent (if you want to work on the whole raster, don't modify anything);


6. Choose if using the current resolution (I suggest to leave it as default);


7. Set the max number of columns and rows (in your case, it should be 20000 columns and 2000 rows);


8. Press the OK button.

For example, the using of the parameters in the above dialog on this sample raster (the parameters I set are chosen randomly):

will generate 100 tiles in the path specified at step 4:

Loading them in QGIS, they look like this:

As @bugmenot123 correctly said in the comments, the result looks weird just because the style of each image fits itself to the distribution of values per image (but the data are perfectly fine).

---

Solution using Python (GDAL)

Another way to obtain the same result is the using of GDAL (gdal_translate).

With reference to the same example described above, you may use this script:

import os, gdal

in_path = 'C:/Users/Marco/Desktop/'
input_filename = 'dtm_5.tif'

out_path = 'C:/Users/Marco/Desktop/output_folder/'
output_filename = 'tile_'

tile_size_x = 50
tile_size_y = 70


ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
xsize = band.XSize
ysize = band.YSize

for i in range(0, xsize, tile_size_x):
    for j in range(0, ysize, tile_size_y):
        com_string = "gdal_translate -of GTIFF -srcwin " + str(i)+ ", " + str(j) + ", " + str(tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(i) + "_" + str(j) + ".tif"
        os.system(com_string)

You obviously need to adapt the values to your specific case.

symbology - ArcGIS 9 doesn't recognize ArcGIS 10 .lyr files

I symbolized a parcel shapefile in ArcGIS 10 with about 2 dozen different colors. I need to use this in ArcGIS 9.3 as well, so I saved it as a .lyr file. This works fine in 10. But when I go to 9.3, I can add the shapefile, but not the .lyr file, or import it for symbolization (geometry doesn't match error). Is there a way to make a .lyr file in 10 so that 9 will recognize it? The .lyr file comes from the same shapefile, so I don't know why it wouldn't recognize it. I want to use this symbolization for a half dozen shapefiles a/o feature classes, so really don't want to manually recreate it. Thanks.

intersection - PostGIS > Intersect two circles. Each Circle must be built from Long/Lat degrees plus radius

I'm trying to figure this out, so I wanted to try asking to the experts.

This question is related to PostGIS (from postgresql).

What I want is:

1) Based on a longitude and latitude I want project a single point (similar on how you can pass the lat and lon to a google map and then have a single point projected).

2) Create a circle geometry around that with a radius. (I'm going to determine the radius based on the horizontal accuracy of the Lat/Lon)

3) Then, I want to do 1 and 2 again to another Lat/Lon combination. And see if those two circles touch.

---

I suppose for point 2), I could use a constructor like ST_GeomFromText, however I dont quite understand yet how I can project a single point based on two spherical angles (lat/lon).

Then for point 3, I guess I could use something like ST_Touches

Answer

I'll post this code to solve part of the problem, as I don't know what you're trying to do with the map part, but here goes:

SELECT 
ST_Touches(temp.point1, temp.point2) as geom_touches
, ST_Intersects(temp.point1, temp.point2) as geom_intersect

FROM (

SELECT * FROM 

ST_Buffer(ST_Transform(ST_GeomFromText('POINT(-105.05083 39.74823)', 4326), 2877), 1500) as point1
, ST_Buffer(ST_Transform(ST_GeomFromText('POINT(-105.04428 39.74779)', 4326), 2877), 1500) as point2

 )as temp

This will:
- create 2 geometries in a temp table, both in SRID:4326 using ST_GeomFromText
- transform those two geometries to a projected coordinate system in Feet (state plane Colorado)
- buffer those points by a distance (1500 feet)
- use an outer SQL statement to test if they touch, and if they intersect

And the results look like this:

And the geometries look like this:

arcgis desktop - Why do ArcToolbox tools give ActiveX error?

I work with ArcGIS 9.3. When I try to run any ArcToolbox tool I receive the message:

Window Internet Explorer [Title of the Window]

One or more ActiveX controls could not be displayed because either:

1. Your current security settings prohibit running ActiveX controls on this pages, or



2. You have blocked a publisher of one of the controls

I tried to solve the problem but I could not. Is anybody who has an idea what is going on?

gdal - Create MBTiles from geojson in raster format

I am trying to create raster file in the Mbtiles format.

I tried ,

$ gdal_rasterize -ot Byte  -burn 255 -burn 0 -burn 0 -ts 4950 4090 -of MBTiles station.geojson station.mbtiles 

I'm using,

$gdalinfo --version : GDAL 2.1.3, released 2017/20/01

Error:

> ERROR 1: Cannot find min and max zoom_level
> 
> ERROR 6: Could not find an appropriate zoom level that matches raster
> pixel size
> 
> ERROR 6: Only EPSG:3857 supported on MBTiles dataset
> 0...10...20...30...40...50...60...70...80...90...100 - done.

> 
> ERROR 6: IWriteBlock() not supported if georeferencing not set
> 
> ERROR 6: IWriteBlock() not supported if georeferencing not set
> 
> ERROR 6: IWriteBlock() not supported if georeferencing not set

EDIT:

$ ogrinfo station.geojson -so -al

INFO: Open of 'stations.geojson' using driver 'GeoJSON' successful.

Layer name: OGRGeoJSON

Geometry: Point

Feature Count: 86

Extent: (-77.272622, 38.766522) - (-76.842038, 39.119927)

Layer SRS WKT:

GEOGCS["WGS 84",

DATUM["WGS_1984",

   SPHEROID["WGS 84",6378137,298.257223563,


       AUTHORITY["EPSG","7030"]],

   AUTHORITY["EPSG","6326"]],

PRIMEM["Greenwich",0,

   AUTHORITY["EPSG","8901"]],

UNIT["degree",0.0174532925199433,

   AUTHORITY["EPSG","9122"]],

AUTHORITY["EPSG","4326"]]

name: String (0.0)

marker-color: String (0.0)

marker-symbol: String (0.0)

line: String (0.0)`

Answer

I'm not certain that this is the case, but GeoJSON is not supposed to be any coordinate reference system (CRS) except for WGS (EPSG:4326).

I think you'd better convert it to a Shapefile first. While you're at it, convert it to the right projection:

 ogr2ogr -f "ESRI Shapefile" station.shp station.geojson -s_srs EPSG:4326 -t_srs EPSG:3857

Now what worked for me was rasterizing to TIFF first:

gdal_rasterize -ot Byte -burn 255 -burn 0 -burn 0 -ts 100 100 -of MBTiles station.shp station.mbtiles

(Note I changed the resolution... play with the -ts parameter.)

And now convert that tiff to MBTiles:

gdal_translate station.tiff station.mbtiles -of MBTILES

Output from my dataset of points across Tasmania loaded fine in QGIS, and ended up with this demonic view, strangely:

Need Projection Help to Buffer Lines in QGIS

I need to create buffers (1km, 3km, and 5km) around the rivers of this shapefile from NOAA, something I remember being relatively simple in my GIS class.

However no matter what I do the buffers I create cover the entire United States. I've tried a number of different CRSs (for the layer I did Right-Click>Save as and set the CRS I wanted to try, I also set the project CRS not just the layer), but the buffer always covers most, if not all, of the continent. Additionally if I use the distance tool to measure this shapefile from coast to coast, it comes out in meters. I could be wrong, but I'm fairly sure that the country is larger then a couple dozen meters across.

qgis - Why is the SelectionChanged signal triggered every time i click on the layer?

i'm making some experiments with the SelectionChanged signal of a QgsVectorLayer.

i try this piece of code :

curLayer.selectionChanged.connect(self.test)

and "test" print a warning message.

The probleme i have is that the signal is triggered every time i click on the layer, even when there's no change of selection ...

Is there something i'm missing in this SelectionChanged behaviour ?

Thanks for your help.

Answer

Ok then, after a few searches :

U got to use the two lists provided by the SelectionChanged signal :

QgsFeatureIds selected,QgsFeatureIds deselected (as mentionned in the API doc...)

to check what is really changing in your selection, of course multiple features selection is to take in account ...

Advantages of using spatialite over shapefile?

I found spatialite more useful than shapefiles as it does not have the limitations of shapefile and it is also portable. Many people here use shapefiles to exchange data and even the experts do not know about this new format.

What are the advantages of using spatialite over shapefiles?

Can it be used instead of shapefile?

Please focus only on those formats which are portable, i.e. can be exchanged using USB sticks. GML, GeoJSON, KML, CSV are not an option and they are not directly editable in GIS.

UPDATE: It has been more than 5 years and the new development is directed towards geopackage which is related to spatialite.

So now question is more like Advantages of using GEOPACKAGE over GEODATABASE?

Answer

Shapefiles are the lowest common denominator of GIS vector data file exchange: send an archive of shapefiles, and you can pretty much guarantee that someone will be able to build a basic GIS from it.

SpatiaLite's advantages include:

• everything's in one file; none of the shp/shx/dbf/idx/prj per layer mess.


• logic as well as data can be included, in the form of VIEWs and TRIGGERs.


• built-in spatial indices, which allow rapid searches of large areas.


• they are real (if slightly limited, mostly in multi-user access) database systems, with no database admin skills required.

But there are some disadvantages:

• not everyone can use them.


• they are still mostly limited to geometries built from nodes and straight lines; if you need splines and surfaces, look elsewhere.

carto - Automatically coloring a large number of polygons in CartoDB

I see that I can use the CartoDB wizard to color polygons I have based on a column value. This works great except that it only colors 6 or so polygons and the rest are lumped as “others”.

I see that I can edit the CSS and add the other polygons one at a time. My map has about 250 of these, so doing it one at a time isn’t practical. Is there a more automated way to, say, randomly assign a color from the palette to all the polygons based on a value in the table? It’s fine if there are duplicate colors.

Answer

I wrote a little script for this that solves the problem by generating all the CSS properties that you can then cut and paste into CartDB:

for (rep = 1; rep <= 800; rep++) {
var color = Math.floor(Math.random()*16777215).toString(16);
if (rep < 10)  {
document.write('#mydiv[territory="T00' + rep + '"] {
 polygon-fill:#' +color+';
}
');
}
  if (rep >=10 && rep <100)
 {

document.write('#mydiv[territory="T0' + rep + '"] {
 polygon-fill:#' +color+';
}
');
}
if (rep >= 100)  {
document.write('#mydiv[territory="T' + rep + '"] {
 polygon-fill:#' +color+';
}
');
}
}

The extra if statements were necessary in my case because I had territories formatted like this "T001" and "T099".

coordinate system - Why not measure elevation from the center of the earth?

In trying to understand Geoid vs Ellipsoid height, I ran across this page.

One of the claims is:

Zero elevation as defined by Spain is not the same zero elevation defined by Canada

Basically, they are indicating that MSL is not constant and that makes comparing things harder.

So my question is: why not measure elevation as a radius from the center of the earth?

I guess I'm proposing {lat, lon, radius} to define a point on or near the earth. This could eliminate many problems in my thinking. (Let's use inertial center of the earth as our reference point...)

Who cares that the topography of the earth is not constant? It wasn't constant in the first place, and to get height above ground, you just need some terrain data and a simple subtraction. I guess I'm proposing using one common reference point that doesn't change.

Answer

One thing to keep in mind is that lat/long is geodetic and not geocentric:

If we were to calculate elevation as a radius from the center of the ellipse, our elevation lat/long would be different than our horizontal lat/long!

This is why there are two different datums. The horizontal datum is just a smooth ellipse, because it's easier to do trig functions on. Most (all?) vertical datums use a geoid, which takes local gravity into account. Establishing horizontal locations has always been a separate process than establishing an elevation.

If you really want to combine horizontal and vertical datums, then...

Welcome to the ITRF fan club! This is the datum and projection for you! One difference is that ITRF doesn't use lat/long at all, it goes all-out Cartesian X/Y/Z to point to any place on or near Earth. This is handy, since the GNSS (GPS, GLONASS, etc) satellites rotate around the Earth's center of mass, it makes a nice, native coordinate system for getting locations from satellites.

Connecting to SQL Server Geodatabase using ArcObject in C#

I am going to save a polygon that is created from points to remote geodatabase. I am new with ArcObjects and can't understand how it works... at this point I want to open a connection to database using ArcObjects and save geometry in it. I am using this code:

public IWorkspace open_ArcSDE_Workspace(string server, string instance, string user,
string password, string database, string version)
{
    // Create the workspace factory.
    Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SqlWorkspaceFactory");
    IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance
        (factoryType);

    // Create the connection properties.
    IPropertySet connectionProps = new PropertySetClass();

    connectionProps.SetProperty("dbclient", server);
    connectionProps.SetProperty("serverinstance", instance);
    connectionProps.SetProperty("authentication_mode", "OSA");//The type of authentication to use. Valid values are DBMS and operating systems authentication (OSA). DBMS is the default mode and is not required.
    connectionProps.SetProperty("user", user);
    connectionProps.SetProperty("password", password);

    // Open the workspace.
    //IWorkspace workspace = workspaceFactory.Open(connectionProps, 0);
    return workspaceFactory.Open(connectionProps, 0);
}

But it throws errors:

enter code hereError HRESULT E_FAIL has been returned from a call to a COM component.

Can you help me and give me some examples that I can understand... The only source I have found for ArcObjects is ESRI and I couldn't find ant solutions for it.

Answer

That should work without problems, double check if the parameters that you are passing to the IPropertySet object are right. Are you writing the "serverinstance" property properly? that's the most common problem.

        IPropertySet propertySet = new PropertySetClass();
        propertySet.SetProperty("dbclient", "SQLServer");
        propertySet.SetProperty("serverinstance", "myMachine\\myInstance");

        propertySet.SetProperty("database", "myDatabase"); // Only if it is needed
        propertySet.SetProperty("authentication_mode", "OSA")
        propertySet.SetProperty("user", "myUser");
        propertySet.SetProperty("password", "myPassword");

But also you can use that way:

    Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SqlWorkspaceFactory");
    IWorkspaceFactory2 workspaceFactory2 = (IWorkspaceFactory2)Activator.CreateInstance
(factoryType);


    // Build a connection string.
    String[] connectionProps = 
    {
        "dbclient=SQLServer", "serverinstance=MyMachine\\SqlExpress", 
            "database=MyDatabase", "authentication_mode=OSA" 
    };
    String connString = String.Join(";", connectionProps);

    IWorkspace workspace = workspaceFactory2.OpenFromString(connString, 0);

In the ArcObjects help the user and password are not included in the connection string, but probably you should added it.

But my question is: which ArcGIS version are you using? It seems that in 10.2 the SqlWorkspaceFactory is nomore supported, there is no online documentation yet and in the local help the chapter "Working with sql workspaces" is nowhere to be found.

Searching for the SqlWorkspaceFactory documentation in the local help, the description is subtle: The SDEWorkspaceFactory.Open method should be used to make connections to all databases and enterprise geodatabases.

If you just give a quick look at the documentation, you may missing that it suggest to use the SDE and not Sql one, without any explaination or anything else.

I hope this will help you.

openlayers 2 - How tow use an ArcGISCache layer along with an ArcGIS93Rest layer?

I want to have two layers in my OpenLayers map. The BaseLayer, as well as the overlay layer are being served by ArcGIS Server.

I have a problem with ArcGIS93Rest when using ArcGISCache as basemap.

The problem is ArcGIS93Rest doesn't show up and there is no HTTP request to the service.

Is it possible to use ArcGISCache and ArcGIS93Rest together?

Here is the code :

var baseLayer = new OpenLayers.Layer.ArcGISCache("AGSCache", "http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer", {

    isBaseLayer: true,                     
    resolutions: resolutions,
    tileSize: new OpenLayers.Size(layerInfo.tileInfo.cols, layerInfo.tileInfo.rows),
    tileOrigin: new OpenLayers.LonLat(layerInfo.tileInfo.origin.x, layerInfo.tileInfo.origin.y),
    maxExtent: layerMaxExtent,
    projection: 'EPSG:' + layerInfo.spatialReference.wkid
});
var forestLayer = new OpenLayers.Layer.ArcGIS93Rest("ForestLayer", "http://192.168.1.36:6080/arcgis/rest/services/forest_editing/MapServer/export", {
    layers: "0,1,2",
    TRANSPARENT: true

});

Answer

I was having problems redrawing the Dynamic map service on every pan & zoom. So to solve this problem, I just called layer.redraw(true);

The important code is as follows:

function init(){
    //The max extent for spherical mercator
    var maxExtent = new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34);

    //Max extent from layerInfo above            
    var layerMaxExtent = new OpenLayers.Bounds(

        layerInfo.fullExtent.xmin, 
        layerInfo.fullExtent.ymin, 
        layerInfo.fullExtent.xmax, 
        layerInfo.fullExtent.ymax  
    );

    var resolutions = [];
    for (var i=0; i        resolutions.push(layerInfo.tileInfo.lods[i].resolution);
    }

    map = new OpenLayers.Map('map', {
        maxExtent: maxExtent,
        StartBounds: layerMaxExtent,
        units: (layerInfo.units == "esriFeet") ? 'ft' : 'm',
        resolutions: resolutions,
        tileSize: new OpenLayers.Size(256,256),                
        projection: mercator,
        displayProjection:geographic,
        eventListeners: { //set event handler
                "moveend": moveEvent}

    });

     cacheLayer = new OpenLayers.Layer.ArcGISCache( "AGSCache",
            "http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer", {
                isBaseLayer: true,

                //From layerInfo above                        
                resolutions: resolutions,                        
                tileSize: new OpenLayers.Size(layerInfo.tileInfo.cols, layerInfo.tileInfo.rows),                        
                tileOrigin: new OpenLayers.LonLat(layerInfo.tileInfo.origin.x , layerInfo.tileInfo.origin.y),                        

                maxExtent: layerMaxExtent,                        
                projection: 'EPSG:' + layerInfo.spatialReference.wkid
            });

    layer = new OpenLayers.Layer.ArcGIS93Rest("ArcGIS Server Layer",
    "http://services.arcgisonline.com/ArcGIS/rest/services/Demographics/USA_Median_Home_Value/MapServer/export?transparent=true&", 
            {layers: "show:0,1,2,3,4", FORMAT:"png24", srs:3857},{singleTile: true});

    //add the 2 layers
    map.addLayers([cacheLayer, layer]);


    map.addControl( new OpenLayers.Control.MousePosition() );

    map.setCenter(new OpenLayers.LonLat(-120, 45).transform(geographic, mercator), 4);

    layer.setVisibility(true); //set visibility
}

function moveEvent(event){
    layer.redraw(true); //force refresh

}

A working JS fiddle is available here: http://jsfiddle.net/devdatta/rDfUx/3/

arcgis desktop - Making Map Document (*.mxd) Locked/Uneditable?

I am wondering if its possible to create a locked map document?

I have a document that multiple users will be accessing and as a way of preventing them from making any accidental errors, am wondering if I can create a map document that is locked as to prevent any modifications.

Answer

Right click the mxd in windows explorer, choose properties, check Read-Only.

This should force 'Save-As' when the user goes to save the map leaving your original in-tact. If you need to make changes to your "template" you'll have to disable this checkmark.

arcgis desktop - How to implement bivariate Ripley's K function?

The attached image shows a forest gap with red pine represented as circles and white pine represented as crosses. I am interested in determining if there is a positive or negative association between the two species of pine trees (i.e. whether or not they are growing in the same areas). I am aware of Kcross and Kmulti in the R spatstat package. However, since I have 50 gaps to analyze and am more familiar with programming in python than R, I would like to find an iterative approach using ArcGIS and python. I am also open to R solutions.

How can I implement a bivariate Ripley's K function?

Answer

After much searching in the back corners of ESRI documentation, I've concluded that there is no reasonable way run a bivariate Ripley's K function in Arcpy/ArcGIS. However, I have found a solution using R:

# Calculates an estimate of the cross-type L-function for a multitype point pattern.
library(maptools)
library(spatstat)
library(sp)

# Subset certain areas within a points shapefile.  In this case, features are grouped by gap number

gap = 1

# Read the shapefile
sdata = readShapePoints("C:/temp/GapPoints.shp")  #Read the shapefile
data = sdata[sdata$SITE_ID == gap,]  # segregate only those points in the given cluster

# Get the convex hull of the study area measurements
gapdata = readShapePoints("C:/temp/GapAreaPoints_merged.shp")  #Read the shapefile that is used to estimate the study area boundary
data2 = gapdata[gapdata$FinalGap == gap,]  # segregate only those points in the given cluster
whole = coordinates(data2) # get just the coords, excluding other data

win = convexhull.xy(whole) # Convex hull is used to get the study area boundary
plot(win)

# Converting to PPP
points = coordinates(data) # get just the coords, excluding other data
ppp = as.ppp(points, win) # Convert the points into the spatstat format
ppp = setmarks(ppp, data$SPECIES) # Set the marks to species type YB or EH
summary(ppp) # General info about the created ppp object
plot(ppp) # Visually check the points and bounding area


# Plot the cross type L function
# Note that the red and green lines show the effects of different edge corrections
plot(Lcross(ppp,"EH","YB"))

# Use the Lcross function to test the spatial relationship between YB and EH
L <- envelope(ppp, Lcross, nsim = 999, i = "EH", j = "YB")
plot(L)

python - Finding Pseudo Nodes using ArcGIS for Desktop?

Trying to write a python scrip to find and dissolve pseudo nodes for a large dataset.

Anyone have any code to perform such a task?

For example, I have a road line with several end/start nodes (pseudo nodes) created by the cartographers.

Easy enough to run a topology but I have upwards of 4,000,000 +/- pseudo node to get rid of.

Software: ArcGIS 10.2.2 Python 2.7

ubuntu - QGIS Server Capabilities Response doesn't contain any Layers (Fedora 15 64-bit)

I'm having similar issues running QGIS Mapserver under a fresh install of Fedora 15 (64-bit). The XML response from GetCapabilities will not show the individual layers inside the QGIS project, just the 'project' layer. All the layers show up just fine under QGIS (1.7.1). For some reason QGIS Mapserver isn't able to parse my project file or find the shapefiles the same way as QGIS does. Shapefiles are located in my HOME directory. The QGIS project file is in /var/www/wms as per the sample configuration file included in /etc/httpd/conf.d/gis-mapserver.conf

Shapefiles are readable by everyone, so is the project. By the way, everything works perfectly under a fresh install of Ubuntu 11.04 32-bit (QGIS 1.7.0 and QGIS-Mapserver 1.7.0) but for some reason I'm up against a brick wall when it comes to Fedora 15 (64-bit). httpd log files don't show anything out of the ordinary. I tried moving the shapefiles into /var/www/wms but Apache (httpd) didn't like that at all. Any ideas?

Here's the procedure I used under Fedora 15 (64-bit) to get where I am. Further down is the procedure I used for Ubuntu 11.04 (32-bit) which worked.

$su -c 'yum update'
$su -c 'yum install httpd'
$su -c '/etc/init.d/httpd start'

Tested apache by going to localhost in firefox

got a Fedora test page . . .

$su -c 'yum install qgis qgis-python qgis-mapserver'

run qgis to test (installed qgis 1.7.1 and qgis-mapserver 1.7.1)

Copied over QGIS projects and shapefiles from another Linux box into my HOME directory, made sure they were readable by everyone

Followed instructions in README file located here:

/usr/share/doc/qgis-mapserver-1.7.1/qgis-mapserver-README.fedora

$cd /var/www
$su -c 'mkdir wms'
cd wms

copied a qgis project into the new wms folder

restarted httpd

$su -c '/etc/init.d/httpd restart'

http//localhost/wms/name_of_qgis_project?service=WMS&version=1.3.0&request=GetCapabilities

(The colon : after http above is missing because evidently I'm only allowed a maximum of 2 links in this post)

XML response shows only the 'project' layer and none of the shapefile layers inside the qgis project

---

If it helps anyone, here is the procedure I used to get QGIS Mapserver working on Ubuntu 11.04. I started from a fresh install and didn't compile anything from source. This worked for me.

$sudo apt-get update
$sudo apt-get install nedit
$sudo nedit /etc/apt/sources.list

Added to the bottom of the /etc/apt/sources.list file:

deb     http://qgis.org/debian natty main
deb-src http://qgis.org/debian natty main

Then

$sudo gpg --keyserver keyserver.ubuntu.com --recv 1F9ADD375CA44993
$sudo gpg --export --armor 1F9ADD375CA44993 | sudo apt-key add -
$sudo apt-get update

Got an error message:

W: GPG error: http://qgis.org natty Release: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY C2A22E8244865A03

$sudo apt-get install apache2 qgis libfcgi-dev libapache2-mod-fcgid qgis-mapserver

WARNING: The following packages cannot be authenticated! libqgis1.7.0 python-qgis-common python-qgis qgis-providers-common qgis-providers qgis-common qgis qgis-mapserver qgis-plugin-grass-common qgis-plugin-grass

(I installed these packages without verification)

Tested apache by going to localhost in firefox

It should say "It Works!"

Copied over QGIS projects and shapefiles from another Linux box into my HOME directory, made sure they were readable by everyone

Made a new folder inside /usr/lib/cgi-bin and copied one QGIS project into it, along with wms_metadata.xml and qgis_mapserv.fcgi

restarted apache

$sudo /etc/init.d/apache2 restart

Tested GetCapabilities

http//localhost/cgi-bin/new_folder_goes_here/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

(The colon : after http above is missing because evidently I'm only allowed a maximum of 2 links in this post)

Was able to add the WMS Layer from inside QGIS without any problem. However, I still cannot get it to load using ESRI ArcCatalog or ArcMap.

Answer

Got it working under Fedora 15 by moving the shapfiles out of my home directory, and then disabling SELinux. Not the ideal solution, but I'll take it for now. I can't figure out how to make it work otherwise.

postgis - How to find the nearest point projected on the road network?

I would like to know if there is a function at pgRouting that given a point (lon/lat) can give you the projection of that point to the road network .

Answer

There are linear referencing functions in PostGIS to help you out project the location of a point along a line. For example, if you have a road, and a point of interest (POI) that is along the side of the road, you can:

Here is the SQL to extract the interpolated POI on the road:

WITH data AS (
    SELECT 'LINESTRING (50 40, 40 60, 50 90, 30 140)'::geometry AS road,
           'POINT (60 110)'::geometry AS poi)

SELECT ST_AsText(
    ST_Line_Interpolate_Point(road, ST_Line_Locate_Point(road, poi))) AS projected_poi

FROM data;

Returns POINT(44.4827586206897 103.793103448276), which is close to the POI, but projected on the road.

The tricky things you might run into is to locate the closest road LINESTRING and/or, your road network might be a MULTILINESTRING, which needs to be broken down into LINESTRINGs in order to work with the above. Also, linear referencing systems work best with projected (i.e., non-Lat/Lon) data, particularly if your data are far north/south from the equator.

---

Update An easier function that can use other geometries, including MULTILINESTRINGs is ST_ClosestPoint. For example, the closest POI on a road network is found using:

WITH data AS (
    SELECT 'MULTILINESTRING ((20 30, 40 70, 40 110), 
                             (40 110, 70 160, 80 190), 
                             (40 110, 25 118, 10 140))'::geometry AS road_network,
           'POINT (40 130)'::geometry AS poi)

SELECT ST_AsText(ST_ClosestPoint(road_network, poi)) AS closest_poi


FROM data;

(Also, as a side note, I'm using a CTE or "WITH" query to supply data. Your query only needs to use the "SELECT" part.)

Seeking open source software package for Remote Sensing?

Here are many questions with great answers about open source GIS software.

I am wondering, what is the best open source software package for Remote Sensing? I would like to learn it and to use in my work.

I used to work with IDRISI, and I've heard about Erdas and ENVI, but they all are not free. I am looking for a free and powerful leader, like Qgis for GIS or R for statistics. With powerful classification, segmentation, Fourier, filters, PCA, etc.

Can anyone please advise me a good free RS software? What are the capabilities, user friendly or with command line? Do any comparison matrices exist?

Answer

There are a few good ones around:

All with the bonus of being able to be used though the QGIS interface using the SEXTANTE plugin like so http://blog.orfeo-toolbox.org/uncategorized/otb-inside-sextante-inside-qgis

QGIS 2.4.0: dissolve tool renders artifacts / errors

The dissolve tool (vector>geoprocessing tools>dissolve) renders some artifacts on shapes with high number of polygons. I’m using a shapefile with 3222 polygons and I’m trying to dissolve it using the field “municipio”, but when the dissolve tool ends calculating the new shape has some artifacts

Does anyone had a similar problem? How to work around it?

I’ve also used v.dissolve from GRASS and the dissolve polygon from SAGA, but the problem persists. I’ve used the same shapefile on ArcGIS 9.3 with the dissolve tool. None artifact/error was produced.

Edit 1

I'm trying to use the v.clean tool as sugested by WhiteboxDev. However I'm getting this error:

An error has occured while executing Python code:

Traceback (most recent call last): File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\gui\ProcessingToolbox.py", line 181, in executeAlgorithm File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 511, in checkBeforeOpeningParametersDialog File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 361, in checkGrass7IsInstalled File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\ntpath.py", line 96, in join TypeError: object of type 'NoneType' has no len()

Python version: 2.7.4 (default, Apr 6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

QGIS version: 2.4.0-Chugiak Chugiak, 8fdd08a

Python path: ['C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing', 'C:/PROGRA~1/QGISCH~1/apps/qgis/./python', u'C:/Users/CERCOTRICHAS/.qgis2/python', u'C:/Users/CERCOTRICHAS/.qgis2/python/plugins', 'C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins', 'C:\PROGRA~1\QGISCH~1\bin\python27.zip', 'C:\PROGRA~1\QGISCH~1\apps\Python27\DLLs', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\plat-win', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISCH~1\bin', 'C:\PROGRA~1\QGISCH~1\apps\Python27', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\GDAL-1.11.0-py2.7-win32.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win32.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\python_dateutil-2.2-py2.7.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\pytz-2014.2-py2.7.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win32.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\six-1.6.1-py2.7.egg', 'C:\PROGRA~1\QGISCH~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\QGISCH~1\apps\qgis\python\plugins\fTools\tools', 'C:\Program Files\SpiderOak\shell_extension_lib\shared.zip']

Not entirely sure if this error is QGIS related... Any idea why I'm getting this error and how can I solve it?