texts: How can I load GeoJSON from URL to leaflet?

texts:
I want to use leaflet's L.GeoJSON to load the GeoJSON from my URL. I have tested to use leaflet-ajax plugin and it worked.

The codes as follow:

<link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />
<script src="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.js"></script>
<script src="dist/leaflet-ajax/leaflet-ajax.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.2.0/jquery.min.js"></script>

  function onEachFeature_SG(feature, layer) {
      if (feature.properties && feature.PLN_AREA_N) {
          layer.bindPopup(feature.properties.PLN_AREA_N);
      }
  }

  var geojsonLayer = new L.geoJson().addTo(map);
  var popup = new L.popup();

  $.ajax({
  dataType: "json",
  url: "mygeojsonurl",
  success: function(data) {
      $(data.features).each(function(key, data) {
          console.log(data.properties.PLN_AREA_N); //has data

              var centerLat = data.properties.cent_lat;
              var centerLon = data.properties.cent_lon;
              var centLatLon = new L.LatLng(centerLat, centerLon); //this is used to place the popup in the "mouseover" function
              layer.setStyle(default_style);

              if (data.properties && data.properties.popupContent) {
                  layer.bindPopup(data.properties.PLN_AREA_N);
              }

              layer.on('click', function (e) {
                console.log(e);
                e.setContent(data.properties.PLN_AREA_N);
                e.openPopup(popup);
              })

              layer.on('mouseover', function (e) {
                  console.log(e);
                  layer.setStyle(default_style);
                  layer.openPopup(centLatLon);
              });
              layer.on('mouseout', function (e) { console.log(e); layer.setStyle(default_style) });
          }

          geojsonLayer.addData(data);
      });
  }
}).error(function() {});

I want it to be simple by using Leaflet directly.


texts: JoinField_management - error 000339: OIDs missing

texts:
I am trying to join a csv-attribute table to a shapefile based on a variable ("pct16").

This is my code:

path_shp = cd + "ArcGIS/"
path_data = cd + "data/2016/raw/"
arcpy.JoinField_management(path_shp + "voting_precincts.shp", "pct16", path_data + "all_precinct_results.csv", "pct16")

When I run it, I get the following error:

  ERROR 000339: Input C:/.../raw/all_precinct_results.csv does not have
  OIDs

The ESRI help advices to add an OID column to the table. My question is, how do I do that in ArcPy? Or, are there any quicker alternatives for my purpose?

A question that seems to be related is Working around Table does not contain OIDs arcpy Error 000339?. Maybe, it contains useful information to solve my problem but I was not able to judge for myself because I did not understand the question.


texts: is there any online WFS gazeteer client which runs on my web browser or exposes http requests?

texts:
My question is somehow similar to this other thread, but in my case I want to test a WFS service:

Is there a web page where I can test a WMS URL?

I look for an online service which lets me add the url of a WFS gazetter service, and then exposes the layers and features available and lets me make filtered requests to that gazetter (i.e. look for a place name) using a web interface.

I am particularly interested in javascript implementations which do actually make WFS requests from my browser, since the final purpose is to expose the underlying http request sent to the WFS server.

In other words: I am tired of trying to construct the right requests on my own. I am not being able to adapt requests from documentation and examples (there are many differences between WFS versions, namespaces, and so on ...)


texts: I'm attempting to use the timestamp function in QGIS 3.4 to record daily progress on projects my technicians work on in QGIS, something that updates automatically when edits are made in the attribute table. I'm not sure what would work best here to allow me to visualize their daily progress for the purposes of generating project metrics. Let me know what you guys think would work the best.


texts: How to build a route that traverses every edge at least once and begins and ends at the same location?

texts:
I have a collection of lines in PostgreSQL w/PostGIS/pgRouting and have built  topology.

I need to produce a route that starts and ends at the same location, but traverses all of the edges.

Can pgRouting solve this routing problem? And if so, what approach or routing function(s) should be used?


texts: Inconsistent word wrapping in legend of ArcGIS Desktop?

texts:
I have been using word wrapping for my legends in ArcGIS 10.3.  I have had some success with this, but am getting errors in the way the labels are displayed.

Here is a pic of the settings:

The error is in the pdf output.  Sometimes the line will be wrapped, but the next line is missing.

In the following example, Towns and Villages (near top left) is cut off in the pdf output.

Here is the pdf output:

Here is the view from ArcMap:

I have searched online and have not been able to find any information about this error.  Has anyone had this error and discovered a fix?

Here is what I have had success with so far:

Adding a space after the offending line makes it appear in the pdf output
Workaround by adding multiple line descriptions and using those in the legend

I have tried adjusting the gaps in the layout settings.  This does not help.  I have also tried changing the label width to different values.  This also does not help.

I tried removing all other elements from the layout as Polygeo suggested.  Once I removed the other layers from the layout, the legend displayed properly in the pdf output, with the lines broken at the correct place and all text showing.  Then when I added all the layers back to the legend, it continued to display properly.  I need to test this further, but at the moment it appears that turning off all layers except the affected one (towns and villages in my example) and then turning the layers back on will fix the legend.

This is strange behavior.  The first time I opened the map today, I wrapped the labels in the legend and the label was cut off for towns and villages.  Then I turned off all the other layers in the map, and the label was no longer cut off.  I exported the map and the pdf output was fine.  Then I turned all the layers back on and the label was still displaying correctly in the map and pdf.  Then I restarted ArcMap expecting to get the wrapping error again, and I didn't get the error this time.  I'll test this again tomorrow.


texts:
I am trying to find all points within a five-mile radius of a given point. I have a query like this:

SELECT * FROM table WHERE ST_Contains(ST_Buffer(geomFromText('POINT(0 0)', 4326), ?), latlon)

I cannot figure out what I put in place of ? (radius) to get five miles. Everything is in EPSG 4326, and according to PostGIS documentation (as best as I can tell), my radius should be in meters. If I put in 12,070.0m (approximately 5-miles), I get matches half way across the country. Does anyone know what I am missing?


texts:
I used the address range VB script

Function FindLabel ( [TOADDL], [TOADDR] )
  FindLabel = [TOADDL]  & "/" & [TOADDR]
End Function

Tried

Function FindLabel ( [TOADDL], [TOADDR] )
if (([TOADDL] > 0) AND ([TOADDR] > 0)) then:
  FindLabel = [TOADDL]  & "/" & [TOADDR]
End Function

But Still did not work

To display the address ranges but my map comes up with 0s all over the map. What else can I add to the coding that would suppress and not display the 0s on my map. See attached left side.


texts: I am a real beginner with QGIS 2.8. I am able to generate a contour map from a raster DEM, but the contour map is very 'busy' because the land form is pretty flat. There are lots of short contour loops that clutter up the output which would be easy to remove if there was a way to select all contours greater than a certain length. QGIS is able to determine the contour length if I inspect a particular contour line, but this is a 'derived' value, and I do not find a way to use this derived value as a selection criterion.


texts: How do I construct a REST query for ArcGIS to obtain the nearest SSA field office?

texts:
I am trying to use one of the US Social Security Administration's APIs, which is an ArcGIS REST service API. It looks extremely powerful, but overwhelming with how complex it is. I am hoping for an example for a simple query.

I am using Python but might just construct the API using the Requests module rather than learning a new library, so I don't need Python specific info.

Here is the endpoint: http://services6.arcgis.com/zFiipv75rloRP5N4/ArcGIS/rest/services/Office_Points/FeatureServer/1

My goal is to be able to tell an SSA recipient where the nearest 2-3 offices are, based on their address, which can be geocoded to latitude/longitude if necessary (it would be great to use the ArcGIS API for that too. Otherwise, I'm within the free tier of Google's geocoding API). I successfully queried to obtain the office matching a zip code, but it looks like I can provide a Geometry of type Point (lat/lng?) and find offices within a set distance. I'm just not understanding how to provide the information to the API nor where to find this in the documentation.


texts: QGIS Layout/Print Composer: Automatic table refresh upon creation from template

texts:
This Question is related to this one:
QGIS print layout: Pulling the name of the map theme from table
Currently I have two specific templates I use for A3 and A0 plans both of them are structured the kind of the same. But there is kind of a small problem. You see there is attribute table frame that pulls a certain value from a table, the annoying part is to make sure this number shows up I have to refresh the table...manually each time I create a new layout from these templates. I want to have this automaticity upon creation. So how?


texts: GeoServer GWC seeding in 900913 slower than 4326

texts:
Using the integrated GeoWebCahce for GeoServer, I am trying to create a tiled cache for a layergroup and having a problem when trying to seed in EPSG:900913.
Even after increasing both the minimal and maximal heapsize and making sure marlin renderer is enabled seeding for EPSG:900913 only works for me in one process at a time and a very slow process at that.

For example I currently have a cache for levels 1:14 but when i try to seed level 15 (~200,00 tiles for my area), then I can only seed with one process or all my processes end prematurely (like described here.), But if I use EPSG:4326 with no pre-existing cache for all those levels (1:15, ~950,000 tiles) then I can use multiple processes and finish all of them in a few hours.

The EPSG:900913 tiles seed at a rate of ~100-200 every 15 minutes, this is very wierd for me but I need them too.
What else might I be able to do?

I am using an Ubuntu Server 16.04 VM on a windows host with 6GB RAM, GeoServer 2.14 Snapshot.


texts:
I need to generate isochrones considering emergency vehicle’s speed. I have checked several packages in R and different plugins in QGIS. I found QNEAT3 (Iso-Area as polygons) as the only tool that the speed can be configured in but I am not sure if this plugin is the appropriate one in my case or not.
Selecting optimization criteria is mandatory in QNEAT3. It means that I have to choose shortest path or fastest path, while I just need to generate isochrones and it should be created using a feasible path for an emergency vehicle and I am not looking for routing optimization. My question here is that how does QNEAT3 works concerning one of these optimization criteria? Does it consider road network for creating iso-Area polygons or the iso-area polygons which QNEAT3 creates maybe be created through paths which a car could not traverse?
Another question is what is the difference between default speed and speed field?


texts: Calculate average values of points intersecting grid

texts:
I'm using QGIS 2.18.16.
I have a layer with 100k locations, each with several values in the attribute table (ex: mass, E impact, v, ...) for each point (imported via CSV file).
I use a vector grid to count the number of spots for each square of 5m by 5m = ok so far. This gives me a heatmap of the number of hits.

I would like to add to each square an average of another value of each spot. (from csv file and in attribute table of hits layer)
Ex: square A, has 250 hits(= individual values) those are already counted. But those hits in square A have each a different Energy impact. I would like to get the average E impact of these hits in each square for the same vector grid.
This gives me different heatmap in the same vector grid. Not with the number of hits per 5x5m square but with the averag E impact.

The final goal is to generate a final layer by multiplying 2 or more values (number of hits, average E impact, presence of people, ...) for each same 5x5m square to obtain a final kind of heatmap that combines those factors for each 5x5m square.

Can this be done?


texts: Dataset '{}' does not exist or is not supported in arcpy

texts:
I am trying to insert some variable into a simple arcpy call using .format(var1, var2) as I would in a string. I am new to using Python 3.6, so maybe there's been some change in how this is done. When I run the following, I get a

  ERROR 000732: Layer Name or Table View: Dataset '{}' does not exist or is not supported

error, which then causes an Execute Error for bad parameters. What am I doing wrong here?

  select = arcpy.SelectLayerByAttribute_management("'{}'", "NEW_SELECTION", "FIPS_CODE = '{}'").format(name, item)


texts: Are there spreadsheet functions to query route planner for travel time and distance?

texts:
I would like to realise a spreadsheet where the columns are as follows.

Place A
Place B
Distance by road between A and B
Travel time by road between A and B

I thought this might be possible using the 'Google' functions of Google Docs' spreadsheet, but I've found none that does the trick.
In the end, I could knock up an app to do this using the Google Maps API, but I would rather avoid it if I can.


texts: Labels with white (or other color) background

texts:
I tried to give my labels somekind of box and a backgroundcolor (or boxcolor) in QGIS, but i couldnt find the right option.
My labels are not readable because my layers are too dark.
I found offset but would prefere a small box in white.
Can u help me?


texts: QGIS: SAGA-GIS module "Add raster value to point" Error: unknown option 'INTERPOL'

texts:
I encountered an error while using the "Add raster value to point" SAGA-GIS module in QGIS (2.16) to retrieve altitude values from a raster DEM previously interpolated from TIN points (as vector features).

I use SAGA-GIS version 2.3.1.

This bug is fully reproducible for me.

I've selected an other points layer in the first field of the module window.
The raster DEM in the the second field.
The interpolation method was set to "[0] Nearest neighbor".
File saved to /temp
[x] Open output file after running algorithm.

Here is the log:
http://pastebin.com/raw/1F5NGz8b

The same operations are working on SAGA's side.


texts:
I am playing with leaflet plugin active-area (https://github.com/Mappy/Leaflet-active-area), but I am having an issue with autopan, specifically regarding autopan when popups are opened.

Desired Behavior:

When a popup is opened, the map should pan to contain that popup within the active area div.  This is what a leaflet map does automatically, but I want it to apply only to the active area div.  But it's not doing that.  Someone had put an issue in (https://github.com/Mappy/Leaflet-active-area/issues/13), but there's not a lot of detail there, and I'm still having the same issue described.

I have a codeblock which generates random markers, and opens the popup of the most recently generated popup.  Here's my code:

// Make empty array to recieve randomly generated markers
randomMarkerArray = [];

// Useful random number function
function randomNumber(min, max){
   return ( Math.random() * (max - min) ) + min
};

// Define buttons
var randomMarkerButton = document.querySelector('#randomMarker');
var clearMarkersButton = document.querySelector('#clearRandomMarkers');

// Empty variable to recieve the layergroup which is created when a the button is pressed
var randomMarkerGroup;

// When the button is pressed:
randomMarkerButton.addEventListener("click", function(){
   // Get map bounds:
   let mapBoundLeft = leafletMap.getBounds()._southWest.lng;
   let mapBoundRight = leafletMap.getBounds()._northEast.lng;
   let mapBoundTop = leafletMap.getBounds()._northEast.lat;
   let mapBoundBottom = leafletMap.getBounds()._southWest.lat;

   // Make random number within bounds
   let randomLng = randomNumber(mapBoundLeft, mapBoundRight);
   let randomLat = randomNumber(mapBoundTop, mapBoundBottom);

   let position = [randomLat, randomLng];

   // Create a random marker
   let randomMarker = L.marker( position )
  .bindPopup(`This is a randomly placed marker

     Latitude: ${randomLat}

     Longitude: ${randomLng}
`);

   // Give each random marker its own name with its index for easy reference
   randomMarker._name = `Random Marker ${randomMarkerArray.length + 1}`;
   // Push it into the array
   randomMarkerArray.push(randomMarker);

   // Create layerGroup using the array
   randomMarkerGroup = L.layerGroup( randomMarkerArray )
      .addTo(leafletMap);
   // Open popup on most recently added marker
   randomMarkerArray[randomMarkerArray.length - 1].openPopup();

})

// Code to clear the random markers you've generated so far
clearMarkersButton.addEventListener("click", function(){
   // Remove the layergroup of random markers
   randomMarkerGroup.remove();
   // Empty the array to truly refresh the random marker layer
   randomMarkerArray = [];
}, false);

See my codepen:

https://codepen.io/slutske22/pen/mdbGJgq?editors=0010

Behavior I'm seeing

When the first marker is generated, it is properly generated within the active area, but then the map pans away from the marker.  Every marker generated after that is then generated in the active area, but the map does not pan to contain the popup within the active area bounds.  Even more frustrating is when you pan the map, then hit the 'Add Random Marker' button.  The map attempts to pan back to where it was before.  But it is not at all consistent.  I am not sure if my problem is in my code, but I dont think so, as it works just fine without active area.  I put in an issue with the gitHub, but no one responding.  I am somewhat new to coding.


texts: Convert SHP to GPX with predefined SHP style in QGIS

texts:
I had defined graduated style in QGIS for points contained in shp format. I had saved this style in sld and qpj format so now I am able to open shp in QGIS and it automatically loads its style (I am not sure if it would with other GIS software assuming I have sld).

But what I would like to do is to convert/save/export this shp with defined style to gpx, kml or kmz (preferably to gpx).

I have been trying...

Right click on a layer
Save as...
In Save vector layer as...: Format: GPX, Choose File Name, Symbology export: Feature Symbology (and Symbol Layer Symbology), GPX_USE_EXTENSIONS=YES, FORCE_GPX_TRACK=YES + other default options.

...but I havent had any luck.

Can this be done?

My goal is to upload gpx to Google My maps.


texts: Is it possible in ArcGIS 10.x that if I select a feature in one layer it will also select the related feature in another layer?


texts: DXF entities not showing on import into QGIS

texts:
I regularly get sent DWG files containing Ordnance Survey base mapping (variously with & without additional layers). I usually re-scale & reposition these in Draftsight, or just run a conversion in Teigha Converter. So far so good.

I have received a DWG with OS mapping & a site boundary. If I convert that to DXF(I usually have success with 2010 ASCII), all I get in QGIS is the site boundary polyline & a load of points.

If I open in a CAD program (Draftsight,QCAD etc) , the DXF that doesn't work in QGIS, I can see the whole drawing.

There is only 1 polyline entity in the import box, so its not that they're all down at 0,0 or something.

I can reproduce this on Linux & Windows, so not related to that.

Any idea where to look?


texts: Set attributes for newly created features in OpenLayers

texts: I have a vector layer and label: "${title}" is used as one of style property for its rendering. Is it possible to set feature.attribute.title property before feature will be drawn by DrawFeature control? Now for newly created feature I get "undefined" label. Of course I can use "beforefeatureadded" event, but my layer contains many features and I want this event triggers only for newly created features.


texts: How to Serve a Cached/Mixed Raster to Geocortex

texts:
I have never used a Cached/Mixed Raster format before, but I had these air photos given to me and I have to get them served up to our Geocortex App. I've been trying for weeks to covert it to a JPEG or TIFF or anything, but everytime it is crashing ArcGIS or producing a blank raster in GDAL.

Here are the specs on one of the rasters I need to publish:

We have ArcGIS 9.3 server here which cannot publish Cached/Mixed raster types. I've been pushing to get us up to 10, which can publish these raster, but I haven't got the approval yet.

Is there an open source server, that Geocortex can access as a service, that can serve this type of raster?


texts: When I create stream network in ArcGIS, for example, if there is a cell of more than 100 flowing into that, it assumes stream. Does ArcGIS decide to cells more than 100 or 80? If I decide, what are the effective parameters ?


texts: Calculate the maximum minimum distance from a geo object to multipoint?


texts:
I've got a few polygons (they were previously one multipolygon) which seem to insist on having an anticlockwise node order. As such, when running through various operations in SQL, they throw out the following error:

  Microsoft.SqlServer.Types.GLArgumentException: 24205: The specified input does not represent a valid geography instance because it exceeds a single hemisphere. Each geography instance must fit inside a single hemisphere. A common reason for this error is that a polygon has the wrong ring orientation. To create a larger than hemisphere geography instance, upgrade the version of SQL Server and change the database compatibility level to at least 110.

No matter what I try, they seem to remain in an anticlockwise order (which I'm getting by looking at the coordinates when I use the following):

SELECT
    CAST(
        Geom AS VARCHAR(MAX)
    )
FROM dbo.Geometries
WHERE id = 1

I've tried using various conversions and operations, including converting to geography:

SELECT
    GEOGRAPHY::STGeomFromText
        (CAST
            (Geom AS VARCHAR(MAX))
    ,4283)
FROM dbo.Geometries
WHERE id = 1

(which, when interrogated, reveals the nodes in the same order)

And ReorientObject() too, which strangely results in a null:

SELECT
    GEOGRAPHY::STGeomFromText
        (CAST
            (Geom AS VARCHAR(MAX))
    ,4283).ReorientObject()
FROM dbo.Geometries
WHERE id = 1

I've even tried some really weird stuff (like getting the shape's envelope, then the difference between the shape and the envelope, then the difference between the envelope and this new shape), but everything is coming back with the coordinates in the wrong direction. Not to mention that I've also redrawn the shape manually too, in the right direction.

Does anyone have any idea what may be going on, or if there may be a way to reverse the coordinates from an extracted WKT? Shape as WKT below:

  POLYGON ((140.31513731241938 -37.361175421343845, 140.31692577675838 -37.361112938712338, 140.31651306045254 -37.3603397028184, 140.31579569995952 -37.359285256895539, 140.3157269163172 -37.35910563021951, 140.31563847294319 -37.359089989244993, 140.31513731241938 -37.361175421343845))

EDIT - should point out that these are all pretty simple shapes with few nodes - certainly less than 30 per shape. However, I'd love it if there were a reliable method of correcting them, for standard polygons, multipolygons, and polygons/multipolygons with rings.

EDIT - also tried buffering, ST_Union with the first point, and MakeValid.

EDIT - have also tried the solution suggested in the possible duplicate, that still leaves it in the current anticlockwise order.


texts:
I'm looking for some advice on GIS file formats.

At work, we use Trimble GPS to record features in the field so we've been exporting as shapefiles to make maps of sites. I use QGIS and a coworker uses AutoCAD maAutoCADocad map.

But I'm trying to send more map info the other way (making a map with surveyed lease boundaries etc, and sending to field staff either on there Trimble or to online maps like Google my maps.) the process is simple enough using QGIS, but I can't figure out how my coworker using AutoCAD can do the same.

I know AutoCAD can read shapefiles but can it create them?

Is there a better file format that would allow both QGIS and AutoCAD to create and read and write info. Right now he can export his project as a dxf file that QGIS can convert to shapefiles but there is no useful attribute data.

Just wondering if it's something coworker is doing wrong or if AutoCAD is the limitation.

I'm Open to any suggestions, we can buy new software if QGIS isn't the ideal software.


texts:
After I installed QGIS 2.2 via OSGeo4W, as I click the icon of QGIS, then it shows this up,

so I can't open QGIS 2.2

Here's some look like troublesome lines in my installation log,

2014/03/23 13:42:15 io_stream_cygfile: fopen(/etc/setup/setup.rc) failed 2 No such     file or directory
2014/03/23 13:42:15 Current Directory: C:\Users\Heinz\Installation files
2014/03/23 13:42:15 Root Directory: C:\OSGeo4W64
2014/03/23 13:42:15 Could not open service McShield for query, start and stop. McAfee may not be installed, or we don't have access.
2014/03/23 13:42:17 root: C:\OSGeo4W64
2014/03/23 13:42:17 Selected local directory: C:\Users\Heinz\Installation files
2014/03/23 13:42:17 Menu name: OSGeo4W
2014/03/23 13:42:17 net: Direct
2014/03/23 13:42:19 site: http://download.osgeo.org/osgeo4w/
2014/03/23 13:42:19 mbox note: Unable to get setup.ini from <http://download.osgeo.org/osgeo4w/>
2014/03/23 13:43:00 site: http://download.osgeo.org/osgeo4w/
2014/03/23 13:43:00 io_stream_cygfile: fopen(/etc/setup/installed.db) failed 2 No such file or directory
2014/03/23 13:43:01 io_stream_cygfile: fopen(/etc/setup/timestamp) failed 2 No such file or directory

How to fix this problem？

I am working under Windows 8.1 64-bit

Update#1
I follow the suggestion of the answer in the post, How to debug QGIS installation problems on XP?, download Microsoft Visual C++ 2010 Redistributable and then install qgis 2.2 using OSGeo4W package, and here's the same part of the installation log,

2014/03/23 16:17:39 io_stream_cygfile: fopen(/etc/setup/setup.rc) failed 2 No such file or directory
2014/03/23 16:17:39 Current Directory: C:\Users\Heinz\Installation files
2014/03/23 16:17:39 Root Directory: C:\OSGeo4W64
2014/03/23 16:17:39 Could not open service McShield for query, start and stop. McAfee may not be installed, or we don't have access.
2014/03/23 16:17:40 root: C:\OSGeo4W64
2014/03/23 16:17:40 Selected local directory: C:\Users\Heinz\Installation files
2014/03/23 16:17:40 Menu name: OSGeo4W
2014/03/23 16:17:40 net: Direct
2014/03/23 16:17:43 site: http://download.osgeo.org/osgeo4w/
2014/03/23 16:17:43 io_stream_cygfile: fopen(/etc/setup/installed.db) failed 2 No such file or directory
2014/03/23 16:17:43 io_stream_cygfile: fopen(/etc/setup/timestamp) failed 2 No such file or directory

I can successfully launch Qgis browser, but I still have the same dumped message that I can't launch qgis desktop.

Update#2
I have uninstalled the 64-bit version and installed the 32-bit OSGeo4W package. I have checked the installation log, just like the one I post in Update#1, and the rest part of it is fine.

I still got the same error message, but unlike the 64-bit version message, the same mbox would appear after I click the confirm button, this time I just have to click confirm in the mbox one time. And I can still launch qgis browser.

Are there someone successfully launching qgis 2.2 on windows 8.1 x64 computer？


texts: Error with Python Script tool for ArcGIS 10.2

texts:
I am trying to create an ArcGIS map tool that includes creating TINS and uses parameters.  This tool works fine when taking the data from shapefiles, but when I use map layers, which is what I want it stops working at CreateTin tool with the error:
Failed to execute (CreateTin).

Here is the TIN creation code that I am using:

# Import arcpy module
import arcpy, sys, string
from arcpy import env
from sys import exit

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")

# Coordinate System
prjfile = r'C:\Data\project.prj'
Coordinate_System = arcpy.SpatialReference(prjfile)

#Set parameters
Points1 = arcpy.GetParameterAsText(0)
Points2 = arcpy.GetParameterAsText(1)
Polygon = arcpy.GetParameterAsText(2)

# Parameters using elevation points for Create TIN tool
TIN1 = "test1"
TIN2= "test2"
Params1 = Points1 + " elev Mass_Points <None>"
Params2 = Points2 + " elev Mass_Points <None>"

# Create TIN models
arcpy.CreateTin_3d(TIN1, Coordinate_System, Params1, "DELAUNAY")
arcpy.CreateTin_3d(TIN2, Coordinate_System, Params2, "DELAUNAY")

instead of using a workspace I want to do everything in the current map but I am not sure how to use the arcpy.mapping module to read the layers into the map. I was thinking this might be an issue.

Still new to this but I need to learn.


texts: How to filter big datasets in Geoserver WMS?

texts: Is there a good way to filter datasets with approx 50,000 features in Geoserver? I'd have to select between 100 and 200 features based on ids. These ids are different from the FeatureIds. Couldn't find an "IN" operator similar to SQL in OGC Filter Encoding Standard. Combining multiple PropertyIsEqualTo operators performs really poorly.


texts: How to include aggregate field values of clustered points

texts:
I am creating a map that has point features representing things like campgrounds, picnic areas, boat launches, etc.  In cases where many points are clustered in an area (representing amenities for the same park or campground), I'd like to create a single point that includes the values of each point in the cluster, so that I can symbolize it as described in Group and align icons in QGIS atlas.

I created the centroid of the point cluster using ST_ClusterWithin() described in Spatial clustering with PostGIS.

The string_agg function mentioned in How to concatenate strings of a string field in a PostgreSQL 'group by' query seems to fit the bill for aggregating the "amenity" field values.  However, I need to "group by" the spatial cluster of points, not tabular values.

Conceptually, it would look something like below, but the values returned in agg_amenity would be grouped by the points each cluster.  Right now, it just aggregates values from all points.  Concept code is below:

SELECT row_number() OVER () AS gid,
    f.agg_amenity,
    ST_NumGeometries(f.gc) AS numgeom,
    ST_Centroid(f.gc) AS geom
   FROM ( SELECT string_agg(table.amenity, ';') AS agg_amenity,
            unnest(ST_ClusterWithin(table.geom, 250)) AS gc
           FROM table) f;

The input looks something like:

| gid | name   | site_name      | amenity    | geom  |
| 1   | Park A | park camground | campground | point |
| 2   | Park A | park picnic    | picnic     | point |
| 3   | Park A | park swim      | swimming   | point |
| 4   | Park B | park camground | campground | point |
| 5   | Park B | park camground | groupsite  | point |

And the output should look like this:

| gid | name   | amenity                        | geom                        |
| 101 | Park A | 'campground, picnic, swimming' | point (centroid of cluster) |
| 102 | Park B | 'campground, groupsite'        | point (centroid of cluster) |

How can I combine spatial clustering and aggregation of the amenity field (or multiple fields, if possible)?  I am using Postgres 9.5 and PostGIS 2.2.


texts: Making DTM from contour vectors in QGIS?

texts:
I have been using QGIS (2.18.0) for 6 months and now need to make a DEM from contour vectors. From my reading, I know I have to use GRASS.

I am in the process of using this tutorial to achieve that end.

After I run r.surf.contour, I get this:

]1

No error message, but the dtm_1m in not in the layer tree.


texts: What are Differences between Feature and Image Services from ArcGIS Online?

texts:
Esri has recently released what they call their Landscape Layers which look to be mostly image services, rather than map or feature services.

These image services look to have feature level detail where you can change the symbology to many different attributes, like a soil type name for example and even show them in a table.

I believe they are planning to release feature services for each of or most of these layers.

What are the main differences, advantages or disadvantages of these image services vs. feature services?  There looks to be a lot of overlap between the two, or am I missing something here?


texts:
I've been digging around, but I can't seem to find out if it's possible.

I create a series of maps for a client and I need to document all the specifications in the map. So I need to create a spreadsheet/word doc of all the properties set in each layer (i.e. symbology colors, border thickness, definitions queries, label expressions, source, coordinate system...basically all the settings in the layers properties dialogue box) Same thing for the data frame properties settings.

Is there a function in ArcMap to do this?

Thanks,
Mike


texts: QGIS: Join by Location: Option to retain all attributes of spatially related objects of my join layer

texts:
I am facing the problem that I want to join the attributes from one specific column of ALL intersecting features from my join layer to my target layer in case of an intersection (1:n relationship) using QGIS. I can´t find any option to do so having wasted all week on this problem. I am looking for a strategy, concatenating the desired attributes of my join layer in a column of my target layer.

Do you know if there is any way to do so without having to do too much scripting?


texts: Calculate number of people within 50 km of a point

texts:
I have two data sets:

World population (raster) in EPSG:4326, from http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density-future-estimates/data-download
A shapefile with points (distributed around the globe), also in  EPSG:4326

I would like to calculate the number of people within for example 50 km from each point. What would be the best way of doing this in QGIS?


texts:
I have to display an overlay with different colors depending on the zipcode in Europe.
Can I get this information from OSM shapefile ? If yes, how exaclty ? Is there a tutorial somewhere ?

Thanks


texts:
I am an undergrad currently studying GIS. I'm working on a paper and I need real world examples of the ArcGIS Extensions i.e. how have people used them. I need to find articles online but so far haven't had much luck. There are plenty of tutorials though.
I need at least one example article for each extension. I was hoping that someone here could point me in the right direction.

Any links to websites or articles where I could find projects that people have done with extensions?

I need examples of the official ArcGIS extensions. Namely: ArcGIS 3D Analyst, ArcGIS Geostatistical Analyst, ArcGIS Network Analyst, ArcGIS Schematics, ArcGIS Spatial Analyst, and ArcGIS Tracking Analyst.


texts: 404 Encountered Loading Legacy (v2.8) PrintTemplate in ArcGIS JS API

texts:
I've created a simple page which implements the PrintWidget and PrintTemplate functionality

Building this small example using version 3.11 of the ArcGIS API works absolutely fine, however when using version 2.8 (current version I'm working with) PrintTemplate seems to cause a 404 error. Having traced the error it seems unable to find http://serverapi.arcgisonline.com/jsapi/arcgis/2.8/js/esri/tasks/PrintTemplate.xd.js

Running the same file using v3.11 again I can see that it looks for http://js.arcgis.com/3.11/esri/tasks/PrintTemplate.js and this works fine.

It seems at some point from 2.6-3.11 - both the format of the file and it's location have changed.

I have checked and rechecked the JS API for PrintTemplate - it does state that it was implemented in v2.6 so it does not seem like a versioning problem. I have tried to find examples using the legacy/non AMD implementations of PrintTemplate but thus far have been unsuccessful.

Given that I may not be able to upgrade to v3.11 without significant work, is there any solution for this?

Hopefully this is a coherent enough explanation of my problem - again I have only just stepped into the world of ArcGIS so this may prove to be an easy fix.

I'm very new to ArcGIS so forgive me if this is a rudimentary problem.


texts: Determine rasters particular area's within ranges using python

texts:
My directions:

Create a script that determines what areas meet the following
conditions:

Moderate slope — between 5 and 20 degrees
Southerly aspect — between 150 and 270 degrees
Forested — land cover types of 41, 42, or 43

Be sure to use the map algebra expressions of the arcpy.sa module.

I'm getting a final raster image with 3 colors on it when it should only be two...either you are in range of the above constraints or not.

#import arcpy
import arcpy
#set env
from arcpy import env
#set workspace
env.workspace = "D:/UW/Winter_2016/501/w5_more_arcpy8_9_10/ex9/"
from arcpy.sa import *

try:

    if arcpy.CheckExtension("spatial")== "Available":
        arcpy.CheckOutExtension("spatial")
#set elevation raster
        elevraster = arcpy.Raster("elevation")

        sloper = Slope(elevraster)
#Moderate slope — between 5 and 20 degrees
        idealslope = sloper < 20
        idealslope_lower = sloper > 5

        slopeFinal = idealslope - idealslope_lower

        slopeFinal.save("slopeFinal")

#Southerly aspect — between 150 and 270 degrees
        sAspect = Aspect(elevraster)
        s = sAspect > 150
        n = sAspect < 270

        sAspectFinal = n - s

        sAspectFinal.save("sAspectFinal")
#Forested — land cover types of 41, 42, or 43
#forestCover = arcpy.Raster("landcover.tif")

        forestCover = RemapValue([[41,1], [42,2], [43,3]])
        outreclass = Reclassify("D:/UW/Winter_2016/501/w5_more_arcpy8_9_10/ex9/landcover.tif", "VALUE", forestCover, "NODATA")
        outreclass.save("forestReclass")

        properArea = (outreclass & sAspectFinal & slopeFinal)
        properArea.save("properArea")

        arcpy.CheckInExtension("spatial")

finally:
    print "I'm done"

Suggestions? I'm sure there is a much better way to do this but I'm so new just trying to do it step by step, so I can understand the larger thinking process in writing python.


texts: How do I show my WMS layer using QGIS server and OpenLayers?

texts:
I've installed QGIS server and it works fine.
I've put qgs file in the server folder dir and using "Add WMS layer" from QGIS I manage to view the qgs file.

This is how I use GetCapabilities for example:
http://www.localhost.com/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&map=D:/OSGeo4W/apps/qgis/bin/MyLondon.qgs&

Now my problem is that in all OpenLayers WMS examples, adding WMS layer looks like this:

layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",
         "http://vmap0.tiles.osgeo.org/wms/vmap0",
          {layers: 'basic'}
           );

How do I convert my WMS to this format? Where in the URL I set the project name and where the layer name?


texts:
I am looking to attempt to create a web app with the Esri Javascript API that is half map, half scatterplot.  There would be one hosted feature layer in the web map, that has two variables in the attributes.

I would like to mimic some of the functionality that is available with ArcGIS Pro charts:

Scatterplot will respect what is in the current map extent (example)
If you interact with features in the map, they are highlighted in the scatterplot (example)
Vice versa of above (select points in scatterplot, selects corresponding features in the map

Lower Priority Requirements:

Draws a trend line through the scatterplot
Scatterplot will respect the symbology of the hosted feature layer (example)

I have tried building an app with Cedar, but my Javascript skills are very out of date.  Will try to add a JSFiddle once I have something that actually resembles the above.  Would love a headstart from someone who has more experience with Cedar and the Esri Javascript API.

UPDATE

Have made some progress in an application available here.

Remaining Issues:

Selecting a point on the chart = does not select the corresponding feature
Selecting an individual feature in map = does not select point in chart
Getting undefined on the chart points.  Suspect fixing this might fix 1 & 2
Making the chart respect the colors from the classification served up from the hosted feature service not respected  (nice to have)


texts: Generating random networks of a specific area in ArcGIS

texts:
I want to compare the richness of mammals between an implemented biological corridor that connects two patches of forest, against a random generated network connecting these same patches. For this purpose I want to generate a random network of corridors that would connect my two forest patches of interest and but that also has the same area as my "real" biological corridor (that is the one that has been implemented already).

Has someone done this type of analysis in the past and can guide me through on how to do it? Any information on what tools or procedures I could use to generate this type of data will be useful. I am not a spatial analyst expert so I am looking for the easiest way to do this. I have almost all extensions for ArcGIS 10.3.1. except the Network analysist tools.


texts: Examples of geospatial scientific posters

texts: Can you point me to examples of interesting scientific posters featuring results of geospatial analysis and/or geovisualizations?


texts:
The question is the one in the title... is it possible to change the default checkbox in leaflet? in that case, how?
Im using StyledLayerControl, Here is the line of css where i edit the menu with checkbox

    .menu-item-checkbox{
    font-family: 'Open-Sans', Arial, sans-serif;
    font-size: 12px;
  color: #d5e8ee;
}


texts: What coordinates are expressed in 6-7 digits numbers?

texts:
I am processing XML documents containing geographic information, and even though they are coming from the same source they are using different coordinates units. The first one is a regular latitude/longitude:

<linearExtension>
    <linearLineStringExtension>
        <gmlLineString srsName="EPSG:4326">
            <coordinates>12.13538573628832 60.10408278104808</coordinates>
        </gmlLineString>
    </linearLineStringExtension>
</linearExtension>

The second one is different:

<predefinedLocationContainer id="100193" type="PredefinedLocation" version="1">
    <predefinedLocationName>
        <values>
            <value lang="no">Tangenkrysset (55) - Dørdal</value>
        </values>
    </predefinedLocationName>
    <location type="Linear">
        <linearExtension>
            <linearLineStringExtension>
                <gmlLineString srsName="http://www.opengis.net/gml/srs/epsg.xml#32633">
                    <coordinates>180418 6548885, 180291 6548826, 180291 6548826, 180285 6548823, 180218 6548784, 180139 6548707, 180139 6548707, 180078 6548653, 180078 6548653, 179986 6548570
                </coordinates>
                </gmlLineString>
            </linearLineStringExtension>
        </linearExtension>
    </location>
</predefinedLocationContainer>

Are they absolute coordinates? If so, in what system? And if they are relative then in addition to the question of their unit system I should know the reference point which is not present here.


texts:
I have a function to create buffers from everything in a shapefile, but I'm getting an error on some files: Record's geometry type does not match collection schema's geometry type: 'Polygon' != 'LineString'

Here's my function:

def buffer_file(src, dst, buffer):
  if not os.path.isfile(dst):
    with fiona.open(src) as input:
      input_meta = copy.deepcopy(input.meta)
      input_meta['schema']['properties'] = { 'id': 'int'}

      with fiona.open(dst, 'w', **input_meta) as output:
        for i, point in enumerate(input):
          output.write({
            'properties': { 'id': i },
            'geometry': mapping(shape(point['geometry']).buffer(buffer))
          })

What am I missing to make this work with all files? I works with most files I throw at it flawlessly. I just have some problem files.


texts: What are the methods for calculating area from a raster in QGIS?

texts:
I have to calculate area from a georeferenced tiff image. I want to find the whole area in the raster.
The two methods I tried was :
1.Using the area measurement tool in QGIS toolbar.
2.Vectorized the raster . Extracted the number of pixels. Calculated pixel resolution in QGIS. Multiplied the number of pixels with area of one pixel.

The area measured in the 2 methods gave two different values.
Could you tell me any another method other than Calculating area of rasters in QGIS?


texts: How to convert multipolygon to multilinestring

texts:
I try to convert my multipolygon layer to a multilinestring layer. In the postgis codumentation
I found this example:

  (...) -- If you have a table of MULTIPOLYGONs
  -- and want to return a MULTILINESTRING composed of the exterior rings of
  each polygon

  SELECT gid, ST_Collect(ST_ExteriorRing(the_geom)) AS erings
      FROM (SELECT gid, (ST_Dump(the_geom)).geom As the_geom
              FROM sometable) As foo
  GROUP BY gid;

  (...)

My postgresql connection is called database, schema public, default
primary key field gid, default geometry field  geom, my
multipolygon layer poly_layer.

How can I convert the layer poly_layer into a new multilinestring layer
line_layer?


texts: I am trying to copy a polygon feature class from a sde geodatabase and paste it into a file geodatabase but got the message saying "Failed to paste XXX A row with this OID already exists". No same feature class is in the file geodatabase. Also I can copy and paste another line feature class from the same sde gdb to the same file gdb with no problem. Anyone has a clue?


texts: Algorithm for calculating the Central Feature using open source Python

texts:
I want to use an open source programmatic approach through Python for calculating the central feature of a feature class with over 700 points. At the moment the process in the code below works but just very slow. It is just over 2 minutes in comparison to 2 seconds when performed in ArcGIS. The output is the same which is good at least.

The question is; are there any data structures or better algorithms that you can suggest to improve the speed of finding the central feature? I have searched on here and through Google but not returning much information. Has anyone attempted this before? or can anyone put forward an efficient way to find the central feature that trumps what I currently have?

The central feature by definition is the feature whose summed distances to all other features is the shortest.

from osgeo import ogr
from shapely.geometry import Point
from datetime import datetime

start_time = datetime.now()

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\*****\Documents\my_geodata.gdb"
## ope the GDB in write mode (1)
ds = driver.Open(gdb, 1)
## reference the layer using the layers name
lyr = ds.GetLayerByName("my_points")

## massive number that the shortest summed distance sill be less than
shortest__total_distance = 1000000000000.00

## keep track of feature x and y with shortest summed distances
central_x = 0.00
central_y = 0.00

## keep track of the features that have been processed
feature_index = 1

## for each point in the layer
for pnt_from in lyr:
    ## set the total distnace to 0.0
    pnt_total_dist = 0.0
    ## reference the feature at the currnt index
    feature = lyr.GetFeature(feature_index)
    ## access the geometry of that feature
    feature_geom = feature.geometry()
    ## get the x and y coords
    x = feature_geom.GetX()
    y = feature_geom.GetY()
    ## reset the index pointer to the first feature
    lyr.ResetReading()

    ## for each point in the layer
    for pnt_to in lyr:
        ## access the geometry and get the x, y coords
        pt = pnt_to.geometry()
        to_x = pt.GetX()
        to_y = pt.GetY()

        ## calculate the distance between each point and the feature_index point
        pnt_distance = Point(x, y).distance(Point(to_x, to_y))

        ## sum the distances cumulatively for each pair
        pnt_total_dist += pnt_distance

        ## if the total distance goes over the shortest total
        ## then stop for this particular point
        if pnt_total_dist > shortest__total_distance:
            break

    ## if the total distance is shorter than the current shortest total
    ## it becomes the new shortest total and updates the central x and y
    if pnt_total_dist < shortest__total_distance:
        shortest__total_distance = pnt_total_dist
        central_x = x
        central_y = y

    ## reset the index pointer to the first feature
    lyr.ResetReading()
    ## increase the feature index
    feature_index += 1
    print feature_index
    ## when the feature index is > the feature count the process is over
    if feature_index > lyr.GetFeatureCount():
        break

## print the coordinates of the central feature
print central_x, central_y

## This part creates a feature class with the central feature as a point.
"""
## create a new point layer with the same spatial ref as lyr
out_lyr = ds.CreateLayer("central_feature", lyr.GetSpatialRef(), ogr.wkbPoint)

## define and create new fields
x_fld = ogr.FieldDefn("X", ogr.OFTReal)
y_fld = ogr.FieldDefn("Y", ogr.OFTReal)
out_lyr.CreateField(x_fld)
out_lyr.CreateField(y_fld)

## create a new point for the mean center
pnt = ogr.Geometry(ogr.wkbPoint)
pnt.AddPoint(central_x, central_y)

## add the central feature to the new layer
feat_dfn = out_lyr.GetLayerDefn()
feat = ogr.Feature(feat_dfn)
feat.SetGeometry(pnt)
feat.SetField("X", central_x)
feat.SetField("Y", central_y)
out_lyr.CreateFeature(feat)

feat = None
"""

ds = None

print datetime.now() - start_time


texts: Tolerance setting in watershed segmentation(EBimage/lidR) in R

texts:
Can anyone here explain to me what the "tol" parameter in the watershed algorithm implemented in the BioConductor package EBImage (and lidR) does?

I find the explanation in the EBImage help difficult to understand:

  The minimum height of the object in the units of image intensity
  between its highest point (seed) and the point where it contacts
  another object (checked for every contact pixel). If the height is
  smaller than the tolerance, the object will be combined with one of
  its neighbors, which is the highest.

Shouldn't it be: if the height is higher than the tolerance, it will be combined with its neighbor?


texts: Process Landsat8 data in AWS S3 from EC2 instance

texts:
I'm wondering what is the best way of processing huge amounts of Landsat8 images stored in the AWS S3 buckets from an Ec2 instance located in the same availability zone.

Should I download the scenes that I need each time I have to process them and then delete when I'm done, and do the same thing every time I need to do some processing?

Or is there a better way, like mounting the S3 bucket into the EC2 instance? I have seen tools like Fuse for mounting, but I am not sure if this is the best way of processing the data.


texts: Shapefile layers in QGIS don't display unless I zoom out

texts:
I have two shapefile layers in a Qgis project that have problems with zooming. One is a line I copied/pasted in a new layer from a set of height contours and the other one is a polygon layer I created by dissolving all features in another layer.

In both cases, the layers won't display unless the entirety of the layer shows in the canvas. When I zoom in or drag the view, as soon as part of the layer is out of the canvas the entire layer disappears.

In the case of the polygon layer, if I delete it and generate it again it works properly, but when I start filling the rings it reaches a point when it stops working well (I want to end up with a single polygon enclosing all features in the original layer). I guess it must have to do with the geometry of the layer, but I haven't found the cause.

Both layers previously worked well, and removing them/adding them again doesn't fix the issue.

All other layers in the project work well (including the set of height contours and the polygons layer I used to generate the troublesome layers).

This issue has shown in both Qgis 2.4 and 2.6

Is there anything I'm doing wrong?


texts: ArcGIS Attribute Table - Latitude, no Longitude

texts:
I received a shapefile from an online source, and in the attribute table it shows a column and values for latitude; but there is nothing for longitude.

I have not used ArcGIS for some time, and I don't remember encountering this issue in the past.

Furthermore, I am trying to align the above shapefile with one that I created with from a series of UTM coordinates. I have read a bit about the 'Projections and Transformations' toolbox but I have hit a block on choosing the right tool/s to use for this purpose.


texts:
I work with Postgis. I have table

CREATE TABLE centerline
(
  gid serial NOT NULL,
  l_from integer,
  l_to integer,
  r_add_pre character varying(15),
  r_to integer ...
}

Export this table with PGSQL2SHP command

command is

pgsql2shp  -f centerline.shp -h localhost -u redsky ....

and when I imported it  back with command

shp2pgsql -i  -s 4326 centerline.shp

it generates table

 CREATE TABLE centerline
    (
      gid serial NOT NULL,
      l_from double precision,
      l_to double precision,
      r_add_pre character varying(15),
      r_to double precision ...
    }

I do not want double type.  Can Anybody help me?


texts:
I'm using QGIS 2.6.1 Brighton, trying to load the following WMS layer:
http://services.nconemap.com/arcgis/services/Imagery/Orthoimagery_Latest/ImageServer/WMSServer

I can successfully load the WMS in ArcGIS, and have spoken with the service providers, who have recently upgraded the service to work with the latest standard for WMS. I was advised to make sure the DPI mode is turned off, but it seems to have no effect.

When I try to add a new WMS layer, the layer appears in the browser window, when I expand the layer, it processes for a few seconds and then stops. Nothing displays in the expanded view and the layer will not add.

What do I need to do to get this layer to load? Does anyone else have success loading this layer in QGIS 2.6.1?


texts: Using rasterio and rio_color command line problem/question

texts:
I'm an amateur when it comes to the command line.

My problem is that the script I'm trying to run (rio) is in the Python\Python36-32\Scripts folder on my C:drive and I'm trying to run it using images found on my f:drive raid as input, but I cant seem to get the command to recognize the path to the imagery even when typing out the full path F:\imagery\input.tif. This command works fine when input.tif and out.tif are in the same dir as the script (Python\Python36-32\Scripts) but I want to be able to run this script on any directory I choose.

How can I achieve this?

I can run RIO from anywhere as I added the scripts folder to the path, so the problem is getting the command to recognize a different directory/drive for the input output. I already tried setting the F:\imagery directory to the path but that didn't work and my google searches are not turning up much for me as I am not to sure what to search for.

Here is the command:

rio color -d uint8 -j 1 input.tif out.tif sigmoidal rgb 2 0.35

and here is a link to rio (im using command line interface) https://github.com/mapbox/rio-color

and rasterio: https://github.com/mapbox/rasterio/blob/master/docs/cli.rst


texts: Multiple Spatial Joins with GeoPandas in Python

texts:
Is it somehow possible to make several spatial joins with GeoPandas without always saving the outcome and loading it again?

I have to join like six different Esri Shapefiles and it takes really long to save them every time.

My code:

import geopandas

#first join
grid = geopandas.GeoDataFrame.from_file("inputpath")
dem = geopandas.GeoDataFrame.from_file("inputpath2")
sjoin = geopandas.sjoin(grid, dem, how = "left")
sjoin.to_file("outputpath", driver='ESRI Shapefile')

#second join with outcome of first join
griddem = geopandas.GeoDataFrame.from_file("outputpath")
slope = geopandas.GeoDataFrame.from_file("inputpath3")
sjoin2 = geopandas.sjoin(griddem, slope, how = "left")
sjoin2.to_file("outputpath2", driver='ESRI Shapefile')

I already tried to directly use the outcome of the first join for the second join like I did in the following code:

import geopandas

#first join
grid = geopandas.GeoDataFrame.from_file("inputpath")
dem = geopandas.GeoDataFrame.from_file("inputpath2")
sjoin = geopandas.sjoin(grid, dem, how = "left")

#second join with outcome of first join
slope = geopandas.GeoDataFrame.from_file("inputpath3")
sjoin2 = geopandas.sjoin(sjoin, slope, how = "left")

But I got following error code:

File "/anaconda2/envs/IP3/lib/python2.7/site-packages/geopandas/tools/sjoin.py", line 51, in sjoin
    " joined".format(index_left, index_right))
ValueError: 'index_left' and 'index_right' cannot be names in the frames being joined

Is there any solution to make this work efficient?


texts: Classification erosion - semi-automatic classification plugin

texts:
I get an error when working Classification Erosion tool in Semi-automatic Classification Plugin (SCP). I use QGIS v3.2.2 and SCP v6.2.4.

The tool is not changing the pixels of the raster layer and creates this short report:

**An error has occurred while executing Python code:**
*NameError: name 'd' is not defined
Traceback (most recent call last):
  File "C:/Users/Sergio/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\maininterface\erosionTab.py", line 65, in erosionClassificationAction
    self.erosionClassification()
  File "C:/Users/Sergio/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\maininterface\erosionTab.py", line 81, in erosionClassification
    sN = cfg.osSCP.path.basename(str(d))
NameError: name 'd' is not defined*

Classification before erosion:


texts: Join on attributes and write vector or create feature from Oracle SDO with pyqgis

texts:
QGIS 2.14 & 2.18

Step 1 :
I do an extraction with QgsDataSourceURI()

#get geometries
uri = QgsDataSourceURI()
uri.setConnection("svdbocsva003.test.blavla.lan", "1521", "tata", "toto", "pwd")
uri.setSrid ('31370');
uri.setWkbType(QGis.WKBMutiPolygon)
uri.setDataSource("OSIRIS", "OSI_SITES" , "SHAPE",b1query ,"OBJECTID")
uri.uri()
sdoLayer = QgsVectorLayer(uri.uri(), "requete", 'oracle')
sdoLayer.setCrs(QgsCoordinateReferenceSystem(31370))
QgsMapLayerRegistry.instance().addMapLayer(sdoLayer)
sdoFeatures = sdoLayer.getFeatures()

All works fine until this step. The layer is added to the project correctly. Features are displayed.

Step 2 : I try to write this layer to hard disk

error = QgsVectorFileWriter.writeAsVectorFormat(sdoLayer, "C:\\user\\test\\test0.shp", "test",None, "ESRI Shapefile")
if error == QgsVectorFileWriter.NoError:
    print "success!"
else :
    print "error"
del error

I get an error message. The file is written but empty as seen as an error has occurred.

Step 3 : I try to create new vector layer with joined data from oracle attributes table and the geom layer with a common field

writer = QgsVectorFileWriter("C:\\user\\test\\test1.shp","test", fields, QGis.WKBMultiPolygon, QgsCoordinateReferenceSystem(31370), "ESRI Shapefile")

# add a feature
for attr in attributes :
    print(attr[0])
    exp = QgsExpression('ID_WS = {}' .format(attr[0]))
    request = QgsFeatureRequest(exp)
    request.setLimit(1)
    feat = QgsFeature()
    feat.setGeometry(sdoLayer.getFeatures(request).next().geometry())
    #feat.setAttributes([1, "text"])
    writer.addFeature(feat)
# delete the writer to flush features to disk
del writer

QGIS is crashing randomly when i iterate on attributes. OGR Logs gives me Source de données invalide () (invalid data source) in the OGR logs. It looks like the problem. Sometimes it does not crash so i can note the file is empty.

EDIT1 : With a workaround by activating the debugger as explained in this post, QGIS is crashing less frequently.When it's not crashing, i get a feature without geometry with a field FID. By the way, i had defined an output layer processing variable. It was generating error in the log. I deleted it as it was not used.
So now the only error message i get is :

Traitement(1): Cannot find variable: svdbocsva003
Traitement(1): Cannot find variable: svdbocsva003

Where svdbocsva003  is the beginning of the host address. I don't understandwhat what is wrong.

EDIT2 : STEP2 solved

The crash occurs when i want to write a file for a second time. So the overwrite is the cause. No need to be loaded in a project for crashing
The right usage of QgsVectorFileWriter.writeAsVectorFormat is with the overrideGeometryType as  QgsWKBTypes.MultiPolygon in my case and forceMulti as True.
With that it's generating well the vector layer

Step 3 solved in the selected answer, my answer


texts: When should you NOT use a spatial index?

texts:
I am asking this because I was mainly working with Oracle but for the past year I've been doubling with PostGIS and SQLServer 2008. Most spatial functions in Oracle won't work without a spatial index returning the ORA-13226 error:

  13226, 00000, "interface not supported without a spatial index"
  //*Cause: The geometry table does not have a spatial index.  //
  *Action:Verify that the geometry table referenced in the spatial
  operator has a spatial index on it.

To me this makes sense. You run a spatial query=you must have a spatial index. But as far as I understand, neither PostGIS not SQL Serve require this. PostGIS even seems to have functions (_* e.g. _STContains) that EXPLICITLY won't use the spatial index.

So the question is- are there any cases where you should NOT use a spatial index?. Not necessarily whether is a 'take it or leave it' approach i.e. it won't make any difference, but where NOT using the spatial index will impove performance? To me, the last sentence is a contradiction in terms but otherwise why PostGIS would provide these functions?


texts: QGIS 3.10 - Field calculator expression running in Attribute Table and not in Processing Modeler

texts:
I use the Field Calculator in the Attributes table and run the operation : "FieldA" - "FieldB", where "FieldA" and "FieldB" are height distances in meter with a double type (i.e. 2,45). The operation run smoothly.
Now, when I do the same operation using the Field Calculator within the Processing Modeler, it does not run the operation without failure. The message error on the Output Preview in the Expression Editor is reporting Cannot convert "FieldA" to double. So it seems that the coding of the expression is not equivalent in both Field Calculator.
In the Field Calculator within the Processing Modeler, I tried :

"FieldA" - "FieldB"
manual input:fields then @FieldA - @FieldB
to_real("FieldA") - to_real("FieldB")
to_real(replace(to_string("FieldA"), ',', '.')) - to_real(replace(to_string("FieldB"), ',', '.'))

I am surely missing something fundamental in how this is running. Would anyone have encountered the same issue?
Many thanks!


texts: Photorealistic texturing in Building footprints using Ortho rectified satellite image

texts:
Is there any way of photorealistic texturing of 2D Building footprints with height field/3D building/Multipatch using ortho-rectified high resolution satellite image (Building Footprints Created from Same Images).
Basic requirement is to give texture in multiple buildings available in satellite images


texts: Can Facebook help show (statistics) how people move from place of living to place of work?

texts:
Maybe some of you guys know if you or FB can help, or who? Perhaps together with some University?

I got this idea that FB probably holds quite a lot of location information on people (that are using FB with location enabled) during different days of week.

Using this data it should be possible to calculate a (partially representative) statistical view of travel in cities, such as Stockholm, that could help to improve traffic planning, and possibly (which is rarely discussed) how to avoid the need for travel.
The latter however, I guess, would require that people's occupation is identified in some form and mapped to corresponding company locations (incl size?) in various locations.

Regards, Fredde


texts: Low Recall for Radar Image Flood Classification

texts:
I am following the example here to learn how to classify flood images.
I am using a single band of the Sentinel-1 radar imagery. I am able to build and train my model without trouble.

The problem is when I test my model, I have a very low recall for the flood class.

Could someone explain why the recall is so low and how I could improve it?


texts: Running C# code from Python tool for ArcMap?

texts:
It appears that the only way to copy individual rows of data while preserving GlobalIDs in code is by using C# or VB.NET.

I need to run the C# code from the above link within a python tool, and I am curious if anyone knows proper techniques or any examples of doing this. I have come across a couple ways to execute C# code within python, such as this and this, but I wanted to see if anyone could show me any examples or tips on doing this for ArcMap 10.

Also, I know that writing the whole script in C# using arcobjects would probably be the 'best' solution, but my organization would like to keep all of our tools in python.


texts:
I have created a DEM from points using the raster interpolation method, the TIN function, based on elevation of each point.  However, when I try to create a contourline, the results show the contours following the lines/boundaries of the 5m cells I've created using the dem. Am I specifying the right values for cell size, search radius, grid size etc with the interpolation?

Ultimately I want to insert the raster I have created into a DEM based on 2m lidar data, where the dem created using the interpolated grid is to replace for missing values/no data coverage in the lidar data.


texts: QGIS 2.8 topological editing creating vertices

texts:
When using the node tool with topological editing enabled on a group of multiple lines, double clicking sometimes doesn't create new vertices. Instead, the line segment gets selected. When this happens, the only solution is to select a different set of lines or a different tool and then return with the node tool to the set of lines I'm currently trying to edit. The issue seems to be aggravated by larger numbers of lines (it doesn't happen at all when editing a single line).

The sequence I run into is I start to edit a group of lines - the first double-click creates a new vertex but subsequent double-clicks result in the impasse and I have to change focus or tools to reset it.

How can I continue to create new vertices by double-clicking?


texts: QGIS - extract raster values greater than 0

texts:
How can I extract values > 0 from an overall raster?

I have a water level grid and a DTM from which I've derived flood depths (using the raster calculator) but obviously the final result contains negative values where there is no flooding.

I want to remove these values (or set them to null) just leaving the positive flood depths.

Many thanks


texts: How to provide automatic overview map in ArcGIS Desktop?

texts:
This is a slightly badly phrased question I know, but I'm not really sure what to call this. My supervisor recently showed me a printed out map that he had produced with a high-resolution satellite image on one side, and then an overview map of the whole of the UK on the other side, with a dot marking the location of the satellite image.

He said that this was produced automatically using ArcGIs, and he didn't have to manually add a point to the overview map with the right co-ordinates.

Unfortunately he has been unable to find time to help me with this, and I have been looking for how to do this, but as I don't know the name of this process it is very difficult to find it.

Does anyone have any ideas what this process/technique is called, and how I should go about doing it?

I'm using ArcGIS 10.


texts: About getting area from DEM below 5000 and below 3000 m elevation?

texts: I want to get the area from the DEM which is below 5000 and 3000. Actually I need those two areas for my Regression Model to calculate the discharge. In that case, how can I find out using python. I need to calculate in each cell of the stream segment. I have tried with Surface Volume tool but it is taking too long time as i need to calculate for all the cells along the stream segments. Is there any easy and fast way to get this information.


texts: Developing real-time navigation App using QGIS API?

texts:
I want to develop a navigation application based on QGIS using C++, for live tracking. I think it is quite normal for GIS to do such things, especially QGIS as it is a very mature GIS platform. But after I searched these key words on gis.stackexchange and osgeo-org, it seems that nobody has developed such an application.

I'm novice for GIS, So I came here looking for help. The application I want has these features:

Vector images can be symbols for target. e.g. I drive from NY to LA, a vector image(.svg/.eps) stands for me.
Real-time displaying. It is necessary to display with the real-time data for real time positioning.
The track should be shown on the map.

I seek anything about how to do such a thing  e.g. websites, examples etc.


texts: Get distance between point and nearest polygon

texts:
I want to be able to compute the distance between a point and the nearest polygon from a shapefile. To be more specific, I have some point coordinates and I want to find their distance from a worldwide coastline polygon shapefile. And I have to do that with GDAL in python.

Could you point me to some direction? I have seen a lot of answers, but since I am a novice on GDAL I get more confused every time.


texts: Does anyone know if I can use Grass to open a .jp2 image file? I am aiming to do some image classification on the image and when I try to import the image into the program a stream of errors comes up in the side window.


texts: Creating stacked pipeline profiles using water utility network tools?

texts:
I downloaded the water utility network editing tool to create profiles of our sewer and storm systems including manholes with elevations.

The profiling tool is pretty easy to use as you just click the manhole at the ends of the segments of interest.  The question is, can I stack two profiles to create a profile that shows both sanitary and storm pipelines?


texts: problem to execute overlay analysis for multiple layer in python script

texts:
I wrote the code to get a series of weighted overly result for multiple layer in python script. But it showed some error massage. What can be the solution?

import arcpy
input1 = "D:\work\helath.gdb\Reclass_Heal1"
a = arcpy.da.SearchCursor (input1, 'Value')
input2 = "D:\work\helath.gdb\Reclass_Heal2"
b = arcpy.da.SearchCursor (input2, 'Value')
input3 = "D:\work\edu.gdb\Reclass_Educ1"
c = arcpy.da.SearchCursor (input3, 'Value')
input4 = "D:\work\edu.gdb\Reclass_Educ2"
e = arcpy.da.SearchCursor (input4, 'Value')
output = r"C:\Users\d\Documents\ArcGIS\Default1.gdb"
input5 = [input1, input2]
input6 = [input2, input4]
value1 = [a, b]
value2 = [c, e]
for i in range(0,len(input5)):
        for j in range(0,len(input6)):
                for k in value1:
                        for l in value2:
                                weighted = output + "\\didar_"+ str(k)
                                arcpy.gp.WeightedOverlay_sa("('input5[i]' 50 'value1'(value[k] 1;value[k] 2;value[k] 3; value[k] 4; value[k] 5; value[k] 6; value[k] 7; value[k] 8; value[k] 9;NODATA NODATA); 'input6[j]' 50 'value2' (value[l] 1; value[l] 2; value[l] 3; value[l] 4; value[l] 5; value[l] 6; value[l] 7; value[l] 8; value[l] 9; NODATA NODATA);1 9 1", weighted)

The error massage is:

  Runtime error  Traceback (most recent call last):   File "",
  line 20, in    File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing_base.py", line
  504, in      return lambda *args: val(*gp_fixargs(args, True))
  ExecuteError: ERROR 999999: Error executing function. Failed to
  execute (WeightedOverlay).

What is the meaning of ERROR 999999 and how this can be solved?


texts: Seeing lat, lot coordinates in raster and getting pixel value

texts:
I am new to QGIS and I have a simple question.

How can I get the pixel value at a given latitude and longitude coordinate in a satellite raster image?

Related to this, I cannot see the latitude and longitude coordinates anywhere. Instead I see at the bottom some x,y coordinates with values like : 411243,2304551.
I am using QGIS 2.01 (last vewrsion) and the image is a Landsat 7 ETM+ in geotiff format.

Many thanks


texts: Converting CSV to DBF using ArcPy giving ERROR 000732

texts:
My code
import arcpy, os, shutil, re
from arcpy import env
arcpy.env.workspace = "R:\GIS_Data\Earthquake_Data\Updates\Eq_updates.mdb"
path = r"R:\GIS_Data\Earthquake_Data\Updates"
In_path = r"R:\GIS_Data\Earthquake_Data\Updates\OGS_Auto_Updates\csv"
Out_path = r"R:\GIS_Data\Earthquake_Data\Updates\OGS_Auto_Updates\dbf_new"
#gdb_name = r"R:\GIS_Data\Earthquake_Data\Updates\Eq_updates.mdb"
connect = "\\"

# convert csv to dbf in a separate folder
# Local variables inputs and outputs: filenames

for filenames in os.listdir(In_path):
    arcpy.TableToTable_conversion(filenames, r"R:\GIS_Data\Earthquake_Data\Updates\Eq_updates.mdb", os.path.splitext(filenames)[0] + ".dbf")
    print("Successfully converted" + os.path.splitext(csv)[0] + ".csv to " + os.path.splitext(csv)[0] + ".dbf")

The error:
File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\conversion.py", line 2249, in TableToTable
raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Rows: Dataset 2000-2009.csv does not exist or is not supported Failed to execute (TableToTable).

Am I missing a import or other command?


texts: How to load GeoJSON files directly onto QGIS?

texts:
Possible Duplicate:
  Plugin or script to get QGIS to display geoJSON in a layer

Is there a way to do this? Or is there a way i can put it onto PostGIS and then load it onto QGIS? I actually used ogr2ogr to load the GeoJSON to Postgis. After that I used ogr2ogr again to pull it out as shapefiles (Using SQL conditions to seperate different geometry types since shapefiles can only hold one geometry type). I then loaded these shapefiles onto QGIS. I want to know if i can directly achieve this without using ogr2ogr or even PostGIS?


texts: Layout toolbar disabled/greyed out in ArcMap?

texts:
For some reason, when I removed and replaced toolbars  in an .mxd document the Layout toolbar is greyed out.

I've tried opening/closing the document and ArcMap, in addition to re-adding the toolbar several times, and also tried going into Data Frame Properties and ensured that automatic radio button was checked.

Any suggestions?


texts: Checking if feature class locked in ArcGIS using ArcPy?

texts:
I'm trying to find a way to determine whether a feature class is locked before trying to alter it.

In one database I want to remove a feature class completely, and in another I want to add extra fields.  I can't do either if they're locked.

I was expecting to find something in arcpy.Describe() but haven't found anything.  Something like

if arcpy.Describe(myFC).isLocked:
    print "Can't proceed - feature class is locked"
else:
    arcpy.Delete_management(myFC)

Is there a way to check feature class locks using arcpy?


texts: Seeking QGIS equivalent to ESRI story maps?

texts:
I am working in QGIS and would like to continue doing so, but I need to create a "story map" that is not time sequenced (as most examples are at Open source equivalent of Esri story maps?)

I did find a 2016 video that references a story plugin- but that no longer exists (https://www.youtube.com/watch?v=9MGupLsU9Eg).

Is there something that I am overlooking?


texts: When to use GiST and when to use SP-GiST index?

texts:
I'm aware of the documentation for both of the indexes, but can't decide which one to use.
PostGIS documentation on these indexes
Example:
Let's say we have table data with the POINT attribute that contains coordinates.
And let's also say we want to preview all points on the map to create a marker clustering. (Only inside of the boundaries of data view). As in the figure below.

This data is not balanced I suppose? what does that mean?
So why to use GiST or SP-GiST?


texts:
I am working with Python with GDAL and Mapfile and I used the following cmd line in order to remove NoData:
cmd = ["gdal_edit.py", src_img , "-unsetnodata"]

And then in my application with my mapfile I defined into my layer:
OFFSITE 0 0 0

In my app, I can see that the NoData is transparent, however I also have "no data holes" into my image.

With Gdalinfo, I have the 4 corners of my image, I need to differentiate the image and the NoData but I don't really know how to that. I thought to add a TIFF behind the image in order to cover the "no data holes" but do you have any better idea?


texts: Why is sometimes geometryType not returned in JSON from ArcGIS Server layer REST API?

texts:
Does someone know why the  geometryType not always get returned from ArcGIS Server layers?

According to the API documentation Layer / Table the geometryType is returned on layers:

//properties specific to layers only
"geometryType" : "<geometryType>",

But on some layers the resulting JSON looks like this:

"geometryType" : "",

The layers I'm testing on is our own and as far as I can see they are set up in the same manner,:capabilities with mapping and web access enabled. Dynamic layers. Some layers return the correct geometry type and some don't.

Update: Clearing the cache See rest admin does not help.


texts: Hide column in attribute table by using PyQGIS

texts:
Is there a way to hide a column in the attribute table by using PyQGIS?

By hiding I do not mean hiding only the edit widget, I would like to hide the column completely.

It is possible to achieve that by right clicking on the column in the attribute table and selecting 'Hide column' or 'Organize columns...' but I did not find a way how to do it by using PyQGIS.


texts:
a) Whats the name of this projection?

Source:http://www.viewsoftheworld.net/?p=5737

b) Is it possible tu use this projection in qgis?


texts:
I'm looking for implementations of an algorithm called "Available Sky". ArcGIS (Spatial Analyst or GRID) is preferred, but solutions in GDAL, SAGA GIS or others are perfectly acceptable.

The description I have is "a method to quantify the influence of terrain on GPS radiocollar performance by creating a variable called "Available Sky" (Rodgers et al. 1997). ... AS is the proportion of sky available to GPS radiocollar through direct line of site in all directions and at all angles without terrain obstructions (disregarding forest cover). ... locations on mountain tops have high AS values ... conversly locations in valley bottoms are low due to mountain ridges on either sides [lateral obstructions]" -- paraphrased from 'GPS radiotelemetry error and bias in mountainous terrain', Robert G. D'Eon, Robert Serrouya, Graham Smith, Christopher O. Kochanny; Wildlife Society Bulletin 2002.

The paper goes on to describe, in outline only, a process of comparing a base elevation model with a coarser "sky" raster set x100m higher than the highest point of the dem. The process is to calculate direct line of sight for each dem point to each Sky point, arriving at an AS value is the proportion of the total number of sky points visible from that location.


texts: Find points from list that satisfy particular rule AND are within a shape/polygon

texts:
I have a list of points which I have plotted using rule-based symbology as 4 categories.

I want to produce (and export) a list of the points that satisfy a certain rule (e.g. "pUXO" = 'Y') AND lie inside of my shapefile.

All I can find information about is counting the number of points in a polygon, but I don't want every point, only the ones that satisfy a certain rule.

If there's a way this could be done using the Python console that would be great.


texts: From classified image to topology in ArcInfo

texts:
I have some coastal 4 band imagery that I have classified to designate land, water, and undetermined. I want to convert this to polygons with topology, and significantly simplify the polygons.

If I convert the raster to vector, then try to apply a topology, then simplify the polygons, Arc Crashes.

(Raster to Polygon>Select By Attribute x3>Export to FeatureDataset as 3 separate classes>Create Topology in ArcCatalog. Started with only Can't-overlap-self rules, but have also tried Can't overlap each other rules, as well as trying one rule at a time. I have tried validating in Arc catalog, with a validation failed error followed by an ArcCatalog Crash. I have tried not validating until in an edit session in ArcMap, which results in a freeze then a crash. I deleted the topology after each error and started a new one.)

If I try to split these polygons, I get the empty output bug in ArcInfo.

(Raster to Polygon>Select By Attribute x3>Export to FeatureDataset as 3 separate classes> Extract/Split tool using a fishnet as the split features> Processes Polygons, completes without error, copies everything from dataset of input feature to Target Workspace, doesn't split polygons. Have seen threads from older versions having this bug.)

If I convert raster to vector, then simplify, I end up with tons of gaps and overlaps left to fix. (Which defeats the purpose of going to a topology first.)

I've got to think this is a fairly common conversion. Is there a specific trick or workflow that I am missing here to accomplish this? Or am I at the mercy of ESRI bugs?

ArcMap / ArcInfo with spatial analysis version 10.0 SP5.

EDIT:Added specific workflows above. Although I've tried a few other variations using mostly the same tools. The usual save and close between steps. Export as a whole before exporting by attribute.


texts: Get links in HTML popup in ArcMap 10.2 to open in default browser instead of IE?

texts:
I have a field with a URL to a PDF. The field is used for hotlinking. When I use the hotlink tool or click on the link in the Identify window, the link opens in my default browser (Chrome), but when I click the link in the HTML popup it opens in Internet Explorer.

Is there a way to get it to open in the default browser instead? Having it open in IE is rather annoying since it does not handle opening PDF files very gracefully.


texts: Writing If-Then Statement in ArcGIS Pro Field Calculator Code Block using Python Parser

texts:
I have two columns in my shape file with text values called GP and Z which contain abbreviations and a third column called UpOrDown which is meant to hold the binary 0-1 generated from the if-then statement.

I am struggling in a few areas that likely are impeding my own research and resolution of this issue.

What is the highlighted area called? I may be using the wrong syntax here but I'm not sure how to find this out because you don't know what you don't know!

Next, I want to use the GP and Z columns only when they are both specific values. For example, if GP is DT and Z is CG I want to output a 1 into the UpOrDown column but if GP is CIC and Z is A I want to output a zero. The categories are extensive and I'm not sure if there is a way to make a library of all the GP and Z pairs that would be 1 or 0 due to my limited python understanding.

Is there any clarity on this that anyone might provide (for a novice)?


texts:
Still to new to Leaflet but with small steps.... I will get somewhere.

I am loading now with jquery.getjson (and this works fine)

// Load the data
jQuery.getJSON("data/lieux_culture_nantes.json", function(data) {
    displayFeatures(data.features, layers, icons);
    var props = ['nom_comple', 'libcategor', 'commune'];
    fuseSearchCtrl.indexFeatures(data.features, props);
});

But for maintenance the user would like to use a CSV file

omnivore.csv('data/lieux_culture_nantes.csv').addTo(map);

But this does not shows the displayfeatures. Can someone give me some direction how to use the CSV option.


texts: PostGIS: How to distinguish overlapping routes (i.e., LineStrings) based on their direction?

texts:
I have a set of routes which are represented by LineStrings in my PostGIS database. These routes represent journeys that people make. I want to come up with routes that (a) overlap, and (b) have the same "direction" (orientation?). I.e., I want to figure out which journeys have the most in common.

To understand what I mean by direction, please see the figure below. The paths A (red) and B (green) start at the hollow markers and end at the solid markers respectively. They overlap to a considerable extent, but they have different directions and hence don't match.

I can find intersections using ST_Intersection(.). Now I need to filter out paths with opposite directions. How can I accomplish this?

EDIT: for now I have an approach based on angle between vectors from start to end points of the paths, but it can clearly fail if the paths intersect but are not "parallel" in their entirety. So I am still looking for answers.

Thanks in advance!


texts:
Is there a way to connect to PostgreSQL in ArcCatalog without SDE? I came across this post from 2011 that suggested using zigGIS but I can’t find a working link to the program.

A little background: At my work we have a small GIS department (3 people) and our current process is to create file geodatabases for each client. The main issue is that we work in a high-transaction environment and we are constantly locking each other when we work on the same project. This is extremely annoying so we're trying to move to a relational database like PostgreSQL.

While it would be nice to switch over to SDE ESRI Canada gave us a quote that is way too much for a staff of three people (something like 40k install and 20k/year maintenance)

Also if anyone uses PostGIS as your main form of data storage, what GIS program do you use?


texts: Do we need to collect ground control points before flying drones or do the drones collect control points during their flight ?


texts: Converting projected coordinates from bbox into long/lat using R?

texts:
I have a DEM file in QGIS and would like to import the the satellite image for the same region into R.

According to QGIS, the CRS used in the DEM file is:

Layer Spatial Reference System
+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs

Layer Extent (layer original source projection)
696584.9723235532874241,167580.0099993922340218 : 696800.0223235532175750,167790.0099993922049180

To download a satelite image I am using get_map. However this only works with long/lat format and I need to convert my projected coordinates into long/lat.

Is there a way how to easily convert the layer extent to long/lat?

I have tried to do the conversion myself using spTransform. Firstly, I have converted the layer extent from QGIS (upper left corner coordinates and lower right corner coordinates) into bbox format (easting, northings).

?bbox
  # first row: eastings (x-axis, longitude)
  # second row: northings (y-axis, latitude)
eastings <- c(696584.9723235532874241, 696800.0223235532175750)
northings <- c (167580.0099993922340218, 167790.0099993922049180)
coords_bbox <- cbind(eastings, northings)
coords_bbox

> coords_bbox_pts
class       : SpatialPoints
features    : 2
extent      : 696585, 696800, 167580, 167790  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs
>

According to my knowledge, the spTransform function works only on Spatial* and therefore I have converted bbox to SpatialPoints and then used the SPatialPoints as input for the spTransform function.

?SpatialPoints
coords_bbox_pts <- SpatialPoints(coords = coords_bbox, proj4string=CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"))
coords_bbox_pts

library(rgdal)
str(trjct_pts)
?spTransform
longlatcoor <- spTransform(coords_bbox_pts,CRS("+proj=longlat +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"))
coordinates(longlatcoor)
longlatcoor

> longlatcoor
class       : SpatialPoints
features    : 2
extent      : -13.40025, -13.3986, 35.41811, 35.42041  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs

I was expecting to obtain a satellite image of a part of the Swiss Alps, but I have clearly done something wrong, because the output I have obtained is not what I expected.

library(ggmap)
library(ggplot2)

# Set a range
lat <- c(-13.40025, -13.3986)
lon <- c(35.41811, 35.42041)

# Get a map
map <- get_map(location = c(lon = mean(lon), lat = mean(lat)), zoom = 14,
               maptype = "satellite", source = "google")
ggmap(map)

What is the problem?

I am not a geographer, there may be quite likely some trivial issue with my thinking.


texts: Extracting polyline of border between differing polygons using R or GRASS?

texts:
I have a shapefile with polygons with an associated variable that can take different values. I would like to extract the borders between polygons of differing value (red and blue) as a polyline (dark black line).

I would prefer solutions in R but will accept other open source solutions.


texts: What is the procedure to load LAS files in QGIS 2.0.1?

texts:
What is the present procedure to load LAS files in QGIS 2.0.1?

After some searching I found a paper outlining the procedure; however, it was for 1.8.something, and related posts here are all rather dated.


texts: I want to plot mean annual precipitation data for a region. I downloaded the precipitation dataset for 1900-2014 from University of Delaware site in netCDF format. Whenever I try to open the dataset in ArcGIS, it gives me more 1000 (114 years * 12 months = 1368 rasters) of rasters with each raster containing the precipitation data of a particular month. I wanted to know that how can I get an average raster data out of the 1368 rasters.


texts: Uniquely defining shortened address for geocoding

texts:
I've searched for existing concepts to define shortened addresses, but didn't find anything yet that fits my needs.

For an example, I have to define object locations for say desks, where desks are numbered D1..Dn for each floor within a building. I need to geocode and reverse-geocode all such objects.

I could of course carry the whole address for all objects, including street, city, etc. But this is a quite verbose. I'd rather use a shortened unique address identifier.

What concepts of shortened address identifiers do you know, and what are their advantages and disadvantages?


texts:
I'm trying to get the sum of "counts" in a raster. I have the code below which can successfully read the counts of each row in raster attribute and stores the counts in a list provided below.

  [763461] [407207] [232269] [137428] [83909] [55074] [33340] [24168]
  [16747] [11488] [7010] [4323] [4021] [2552] [1497] [775] [72]

However, I want to get the sum of the counts. I am providing the code and error below.

code:

import arcpy
from arcpy import env
from arcpy.sa import *

#To overwrite output
arcpy.env.overwriteOutput = True

#Set environment settings
env.workspace = "C:/Subhasis/Test/7-17-14/test"

#set local variables-STI extraction
inRaster ="sti"
inSQLClause = "VALUE >= 9"

#checkout ArcGIS spatial analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute ExtractByAttributes
attExtract = ExtractByAttributes(inRaster, inSQLClause)

# Save the output
attExtract.save("C:/Subhasis/Test/7-17-14/test/STI-9")

inraster ="sti-9"
rows = arcpy.SearchCursor(inraster,"","","Count","")

for row in rows:
    count = row.getValue("Count")
    s1=[]
    s1.append(count)
    print sum(range(s1))

Error:

Traceback (most recent call last):
  File "C:/Subhasis/Test/7-17-14/test/test.py", line 31, in <module>
    print sum(range(s1))
TypeError: range() integer end argument expected, got list.


texts: How can you define distance in QGIS Multiple ring buffer?

texts: I am trying to create multiple buffers around point data. I need to crate buffers of 10 miles, but I cannot find a way to define this. I have tried using both the multiring buffer plugin and the multi-distance buffer plugin, but both of these are not producing the doughnut style rings, but instead just fill the map. I have limited the attribute table to the 10 ring and this still fills the table, so I can only assume the problem is with not being able to state I want miles. Is there a way to define the distances needed?


texts: GeoServer Point Labels not rendering until Scale = 1 : 545K

texts:
I have a series of points that I have reading an attribute to display a label.
I don't see the labels on my map until Scale = 1 : 1M when SOME of them start to show up. But it's not until I hit Scale = 1 : 545K that they start showing up completely.
Is there something I'm missing in my SLD? I've tried numerous things and cannot get these labels to show prior to these zoom levels.

<StyledLayerDescriptor version="1.0.0"
    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">
  <NamedLayer>
    <Name>Simple point</Name>
    <UserStyle>
      <Title>Points with Labels</Title>
      <FeatureTypeStyle>
        <Rule>
          <MaxScaleDenominator>10000000</MaxScaleDenominator>

          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>circle</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#2234cb</CssParameter>
                </Fill>
              </Mark>
              <Size>10</Size>
            </Graphic>
          </PointSymbolizer>

        </Rule>
        <Rule>
          <TextSymbolizer>
            <Label>
                <ogc:PropertyName>ATTRIBUTE_NAME</ogc:PropertyName>
            </Label>
            <Font>
                <CssParameter name="font-family">Arial</CssParameter>
                <CssParameter name="font-size">12</CssParameter>
                <CssParameter name="font-style">normal</CssParameter>
                <CssParameter name="font-weight">bold</CssParameter>
            </Font>

            <LabelPlacement>
                                  <MinScaleDenominator>0</MinScaleDenominator>

              <PointPlacement>
                <AnchorPoint>
                  <AnchorPointX>
                    <ogc:Literal>0.5</ogc:Literal>
                  </AnchorPointX>
                  <AnchorPointY>
                    <ogc:Literal>0.0</ogc:Literal>
                  </AnchorPointY>
                </AnchorPoint>
                <Displacement>
                  <DisplacementX>
                    <ogc:Literal>0</ogc:Literal>
                  </DisplacementX>
                  <DisplacementY>
                    <ogc:Literal>5</ogc:Literal>
                  </DisplacementY>
                </Displacement>
                <Rotation>
                  <ogc:Literal>0.0</ogc:Literal>
                </Rotation>
              </PointPlacement>
            </LabelPlacement>
            <Halo>
              <Radius>1</Radius>
              <Fill>
                <CssParameter name="fill">#fffff0</CssParameter>
              </Fill>
            </Halo>

             <VendorOption name="conflictResolution">false</VendorOption>
             <VendorOption name="goodnessOfFit">1.0</VendorOption>
            <VendorOption name="autoWrap">75</VendorOption>

          </TextSymbolizer>
        </Rule>
         </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts:
Can we find last added layer in pyqgis? I use Data Source Manager | Delimited Text at the end of the process it creates a layer and adds it layers panel:

However, I can not know what the Layer Name is and I have to know that what it is.


texts: Making Operational Layer from WebTiledLayer and MapImage?

texts:
I'm attempting to use the CMV and put a WebTiled layer in as an operational layer.  I've been able to get the layer into the basemap without issue, but can't get sort it out for the Operational Layer.  I'm trying the same thing for a MapImage, it works in a basic viewer and not with the syntax in viewer.js.  Any suggestions on where to begin would be excellent.

    //Working code as a basemap
mapboxPirates: {
        title: 'Pirates (mapbox)',
        basemap: new Basemap({
            id: 'mapboxPirates',
            layers: [new BasemapLayer({
                url: "http://${subDomain}.tiles.mapbox.com/v3/aj.Sketchy2/${level}/${col}/${row}.png",
                copyright: "mapbox, 2015",
                id: "mapboxPirates",
                subDomains: ["a","b","c","d"],
                type:"WebTiledLayer"
                })]
        })
}

    //An operational layer mapimage that doesn't work
    {
  type: 'mapimage',
  title: 'radar',
  id: 'radar',
  options:{
      extent: { xmin: -8864908, ymin: 3885443, xmax: -8762763, ymax: 3976997, spatialReference: { wkid: 102100 }},
      href: 'http://localhost/goes_ecan_1070_100.jpg',
    visible: true,
  }
}

    //Operational Web tiled that doesn't work
    {
  type: 'tiled',
  title: 'The Web Tiled',
  id: 'theTiled',
  urlTemplate: 'http://{subDomain}.mqcdn.com/tiles/1.0.0/vx/map/{level}/{col}/{row}.jpg',
  options: {
    copyright: "MapQuest © 2012",
    subDomains: ["mtile01", "mtile02", "mtile03", "mtile04"],
  }
}