texts:
Some of the geometry fields that are returned from my PostGIS queries are truncated. When I first found this error any geometry fields with more than a few thousand characters would not show up in the pgAdmin Output Pane. So I changed the pgAdmin >> Options >> Query Editor >> Max Characters per column to 32,500 and now the geometry field shows up but it is truncated and ends with "(...)" like the following example;

MULTIPOLYGON(((
    -79.4660582160837 42.211909028322,
    <<..coordinates>>..
    -79.4673861587629 42.212436007559,
    -79.4673871980406 42.212436150477
   (...)

Any query that ends with the "(...)" fails when I attempt to serialize the result. This occurs on any geometries that have more than a couple thousand characters in them.

I'm pretty sure the issue has something to do with these large geometries because I can successfully run geometries that have fewer vertexes.

When I first found this bug I thought it might be restricted to PostgreSQL 9.2 so I installed version 9.3 and have found the same issue. I also reindexed my table and as far as I can tell the table is not corrupt.

I am also using PostGIS 2.2 bundle for PostgreSQL x64 9.3.


texts: Mapserver - Image Dimensions Don't Match WMS Parameters

texts:
Here is the .map file I am using for a simple map with one layer showing the blue marble world:

MAP
  NAME   "world"
  CONFIG "MS_ERRORFILE" "/home/balajeerc/Projects/gis_project_Build_Dir/gis_project/data/logs/mapservice_errors.log"
  CONFIG "PROJ_LIB" "/home/balajeerc/Projects/gis_project_Build_Dir/gis_project/deps/proj.4-master/share/proj"
  CONFIG "ON_MISSING_DATA" "LOG"
  DEBUG 5
  EXTENT -180.000 -90.000 180.000 90.000
  SHAPEPATH "/home/balajeerc/Projects/gis_project_Build_Dir/gis_project/data/service_data/data/wms"
  IMAGECOLOR 255 255 255
  SIZE 1350 675

  WEB
    METADATA
      "wms_title" "Georbis WMS Server"
      "wms_author" "MyCompany"
      "wms_onlineresource" "http://10.0.0.8:19090/rasters/world/wms?"
      "wms_enable_request" "*"
      "wms_srs" "epsg:4326 epsg:3857"
      "wms_feature_info_mime_type" "text/html"
      "wms_format" "image/jpg"
    END
  END

  PROJECTION
    "init=epsg:4326"
  END

  LAYER
    NAME world
    TYPE RASTER
    DATA "world.tif"
    EXTENT -180.000 -90.000 180.000 90.000
    METADATA
      "wms_title" "BlueWorld Map"
      "wms_srs" "epsg:4326 epsg:3857"
      "wms_server_version" "1.1.1"
    END
    PROJECTION
      "init=epsg:4326"
    END
  END
END

Here is the WMS request I am making to the server:

http://localhost:19090/rasters/world/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fjpeg&TRANSPARENT=false&LAYERS=world&mode=map&CRS=EPSG%3A4326&STYLES=&WIDTH=696&HEIGHT=872&BBOX=303302.12823557947%2C-14401959.121379768%2C8834897.477313813%2C-7592337.145509986

The trouble is that even though I specify a width and height of 696 and 872 respectively, I always get an image with the original image dimensions, i.e. 1350 x 675.

Why is this the case? Like other WMS servers, I am hoping that Mapserver will also scale the base layer into the resolution I am requesting from it.


texts: Integrate Google Street View into Openlayers3

texts: I am wondering if it is possible to integrate Google Street View into my Openlayers3 application. I do not want to show Google Maps, only Street View, since the base maps I am working with are more detailed.


texts:
I am trying to configure JSTS to load using ES6 modules so that I don't have to include the entire library.

There are minimal examples of how to do this online, including using an extend function to pull in functionality, however all attempts at getting the isValid method to work have failed:

import extend from 'jsts/extend';
import Geometry from 'jsts/org/locationtech/jts/geom/Geometry';
import Valid from 'jsts/org/locationtech/jts/operation/valid';
import WKTReader from 'jsts/org/locationtech/jts/io/WKTReader';

    extend(Geometry.prototype, {
        buffer: function () {
            return Valid.isValid(this, ...arguments);
        }
    })

    let reader = new WKTReader();
    var geom = reader.read('POLYGON ((80 300, 280 300, 280 80, 80 80, 80 300), (260 280, 180 200, 100 280, 100 100, 260 100, 260 280))');
    console.log(geom.isValid());

Each time all I am getting back is:
Uncaught TypeError: geom.isValid is not a function

Example JSTS ES6 implementation:
https://github.com/DenisCarriere/jsts-es6-example


texts:
Until now I doing the following steps to get rounded coordinates:

create a new x and y column in the attribute table
use Calculate Geometry
round the column value with the Field Calculator using VBS
round([x],0) and round([y],0).

How do I round the coordinates of a point before insert into attribute table in ArcGIS?


texts: how to send question to QGIS mailing list

texts:
I am not sure but everytime I post any question to QGIS mailing list, I get the below message?
Am i doing something wrong?

You are not allowed to post to this mailing list, and your message has
been automatically rejected.  If you think that your messages are
being rejected in error, contact the mailing list owner at
qgis-community-team-owner@lists.osgeo.org.


texts: Progress Bar gets stuck converting Lines to Polygons using QGIS?

texts:
I'm trying to convert lines to polygon using QGIS 2.8.1.

I use Vector > Geometry Tools > Polygons to Lines.

I always get the progress bar stuck on 90%, it can't get to 100%.

Could someone help me to solve this problem?


texts: Most efficient method to join multiple fields in ArcGIS

texts:
I have as set of shapefiles with approximately 300000 polygons each, as well as corresponding point shapefiles and tables with approximately the same number of row (equal or less rows than the polygons FC).

I need to join about 30 fields from the table of the point FC and from the stand alone table. It works, but it is quite slow.

I've tried different method : AddJoin + claculateField, JoinField and arcpy.da.UpdateCursor. Join Field is the slowest and the Cursor is the fastest but I found it tricky to manage when the two tables are not matching exactly. Nevertheless, I was surprised that a spatial join is sometimes faster than a table join: this makes me think that I a missing something. I saw on a Website that I should work with layer instead of the FC, but I used "makeFeatureLayer" in my script and it didn't help much.

Can I improve this with indexing, with other FC types or with another approach ? Note that this is part of a broader Python script with arcpy, so I need a solution that can be used within this script.


texts: XY coordinate data not lining up with other layers

texts: When I add XY (lat,long) data to the map in the same PCS and GCS as the other layers, the points are added to the map but do not line up with the other layers. I am using GCS: GCS_North_American_1983, and PCS: NAD_1983_UTM_Zone_16N. When I change the XY data to WGS_1984, the points appear on the other layers. However, I am unable to run any measuring tool. I am trying to use the Near tool to measure distance between the XY data and a polygon layer. This does not work, probably because they are in different coordinate systems. I've also tried converted the XY data into a shapefile and that doesn't work either.


texts:
How can I bring into JOSM a Esri Map Service?
For example:

http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer

Googling around I found How can I create an adapter between the ArcGIS REST interface and WMS or TMS? but what if the Esri map service does not have WMS or TMS enabled?


texts: ArcGIS Portal instalation on the Windows OS?

texts:
I want to know if is it possible to install ArcGIS Portal on the Windows OS (other than Linux).

Can ArcGIS Portal be installed without ArcGIS Desktop?


texts:
I'm working on a TDOA application. This is just something in my free time to see if my hardware can accurately do it and it sounded like a cool challenge. Regardless if it looks possible or not, I'm still wanting to try.

So that being said... I wanted to see if someone could help me understand how a location is generated from a 3d set of 4+ receivers when they hear something form an emitter.

In my case, I'm using 4 receivers. I have lat/lng, and alt of all of the receivers. I am logging the times. This is where I start to see variances however in my google research. Some say I need to take a time difference of each receiver. Some say I can calculate the distance. I don't fully understand how it all works. I am find not seeing this with actual numbers, but if anyone could give me a walk through on this. I seek a better understanding of how you get from the timed recordings to a x,y and z positions of the emitter. Everything I find online so far is too vague for me to properly follow along.

And I know what I'm describing isn't trilateration, but seems these two have a lot in common.


texts: OpenLayers 3 - how to redraw a single feature

texts:
I would like to be able to trigger re-rendering of a single feature after a user has clicked on it. (I can catch the click and update the metadata that feeds the styling function, but I am then forced to re-render the entire layer to get the updated styling to show:

kmlSource.dispatchChangeEvent();

Performance is horrible redrawing a an entire large layer when just one feature's style has changed. Is there any way I can redraw just one feature?

Even if the rendering wasn't perfect (for example the redrawn feature was on top) that would be fine.

I note this older question on the previous version of openlayers, but I can't find anything addressing the issue on OpenLayers 3.


texts: WFS feature caching in QGIS hitting feature limit

texts:
I am using QGIS 2 and using a WFS with a 5,000 feature limit.

When I load a dataset without caching the features and it hits the 5,000 feature limit, QGIS seems to give up on reloading data. When I zoom in on the dataset, it continues to request the WFS from the original extent, reloading the data, but not in the new extent.

The upshot of it is that when I load a dataset, hit the limit and zoom in, I still can't see the entire dataset without removing and re-adding the dataset

Is this something wrong with the way I have QGIS set up, or is it something to do with the WFS settings?


texts: Intensity map based on a value of the each record

texts:
If I have geolocalized records with a simple true/false value for a field with title "Do you like flowers?"

How can I draw a map that shows me in which areas people like flowers more than in which areas like flowers less.

If I need to after-cook the data it is not a problem but I don't see how can I draw a kind of a "intensity/heat" map with the areas where people like flowers more. Remember that counting "yes" it not enough, has to be the difference between "yes"es and "no"es for each area.

The table would look like this:

Lat, Long, like flowers
1,   1,    y
3,   2,    n
4,   2,    n
4,   1,    y

And I would like to generate something like this:


texts:
I'm still new to Leaflet and have been trying to add a GeoJSON file to a Leaflet map. I've checked the GeoJSON file using geojson.io and it loads there but not in my Leaflet map. I've looked at the answers here: (Adding local GeoJSON files to Leaflet web page) with no success. I have no idea where I'm going wrong!

Here's the link to my map: http://www.huskey.gisdenver.co/map6.html

1st problem code:

var mydata;
$.getJSON("classes_4326.geojson", function(data) {mydata = data;});
var landcov = L.geoJSON(mydata,
        {style: function(feature) {
          switch (feature.properties.landcover) {
            case 'agriculture': return {color: "#32CD32"};
            case 'grassland': return {color: "#C51C0F"};
            case 'water': return {color: "#1E90FF"};
            case 'wetlands': return {color: "#314004"};
          }
        }, onEachFeature: onEachFeature
}).addTo(map);

As @TomazicM commented, it should have been formatted like this:

$.getJSON("classes_4326.geojson", function(data) {
var landcov = L.geoJSON(data,
        {style: function(feature) {
          switch (feature.properties.landcover) {
            case 'agriculture': return {color: "#32CD32"};
            case 'grassland': return {color: "#C51C0F"};
            case 'water': return {color: "#1E90FF"};
            case 'wetlands': return {color: "#314004"};
          }
        }, onEachFeature: onEachFeature
}).addTo(map);

@TomazicM's sharp eyes spotted that I was missing the closing

});

for the

$(document).ready

function at the end of the code.


texts: QGIS:When exporting a print layout as a .svg, the entire map gets exported

texts:
I have a problem, I'm using QGis 3.2.3, when I make a print layout an put a map on it

, if I use some projections (in this case albert projection focused on canada), then export it as a .svg, the entire map gets exported  resulting in a very heavy vector file that is hardly usable. However if I export it as a picture, it works  and only shows what is inside the frame, furthermore I don't encounter this problem if use some other projection, like platte carree.

What am I doing wrong? Is this a glitch?


texts:
I am trying to run Python code for Grass GIS from this site

Specifically, here is my code:

## Import libraries needed to call GRASS using Python
import grass.script as grass
from grass.script import core as grass
import grass.script as grass
grass.run_command('g.region', raster='elevation')

But I am getting the following error:

---------------------------------------------------------------------------
CalledModuleError                         Traceback (most recent call last)
<ipython-input-13-2ea69f70ebcb> in <module>()
      7 from grass.script import core as grass
      8 import grass.script as grass
----> 9 grass.run_command('g.region', raster='elevation')

C:\Program Files\GRASS GIS 7.5.svn\etc\python\grass\script\core.pyc in run_command(*args, **kwargs)
    419     else:
    420         returncode = ps.wait()
--> 421     return handle_errors(returncode, returncode, args, kwargs)
    422
    423

C:\Program Files\GRASS GIS 7.5.svn\etc\python\grass\script\core.pyc in handle_errors(returncode, result, args, kwargs)
    330         args = make_command(*args, **kwargs)
    331         raise CalledModuleError(module=None, code=repr(args),
--> 332                                 returncode=returncode)
    333
    334 def start_command(prog, flags=b"", overwrite=False, quiet=False,

CalledModuleError: Module run None ['g.region', 'raster=elevation'] ended with error
Process ended with non-zero return code 1. See errors in the (error) output.


texts: GEOTOOLS - FileDataStoreFinder.getDataStore() on jar

texts:
The following code works fine on my Eclipse IDE.

private void getLayersAndDisplay() throws Exception {
        URL imageURL = ImageLab.class.getResource("earthlights.jpg");
        File imageFile = new File(imageURL.toURI());
        URL shapeFileURL = ImageLab.class.getResource("countries.shp");
        File shapeFile = new File(shapeFileURL.toURI());
        URL shapeFileURL2 = ImageLab.class.getResource("Brasil.shp");
        File shapeFile2 = new File(shapeFileURL2.toURI());
        displayLayers(imageFile, shapeFile,shapeFile2);
    }

However, I need to compile this to a jar File. A jar doesn't accept the creation of a new File on its content. The problem is, that in the function displayLayers, the FileDataStoreFinder.getDataStore() only accepts File types to work on. Also, I think this part of the code would give me some troubles too:

AbstractGridFormat format = GridFormatFinder.findFormat( rasterFile );
Hints hints = new Hints();
if (format instanceof GeoTiffFormat) {
    hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
}
reader = format.getReader(rasterFile, hints);

Is there any possible solution to that?


texts: How to translate a pixelCoordinate between different map sizes? OpenLayers 6.3.1

texts:
My starting point for this question is the ability to print a map with a given format&resolution, very similar to the "Export to PDF" example at OpenLayers webpages:
https://openlayers.org/en/latest/examples/export-pdf.html

I am not using the jsPDF library, but using similar approach of drawing to a printCanvas.

My challenge comes from when there is a measurement overlay in the map that should be printed. This is a screenshot from my map:

If I print my map "as is", meaning not changing format or resolution making sure the map size stays the same, everything is fine. I have created a canvas "renderer" of the overlay element, and this works decently.

However, when DPI/resolution or Format is modified, the size of the map changes. And hence I need to reposition my overlay.

My initial thought was that this shouldn't be too difficult, in my mind I wanted to implement it like this:

Get the current position of the overlay in pixel coordinates.
Calculate the ratio for width and height by taking new map size w/h and dividing it on old map size w/h
"Translate" the old pixel coordinates to new pixel coordinates by multiplying them with the ratio
Calculate the new actual coordinates by using the Map -> getCoordinateFromPixel() method
Update the overlay position by using Overlay -> setPosition()

A simple small example on paper:

Old map size: 100w 100h
Old overlay position in pixels coords: 50,50
--- Increase map size by the double ---
New map size: 200w 200h
widthRatio = 2
heightRatio = 2
New overlay position in pixel coords: 100,100

However when I use my translated pixel coordinates, and update the position of my overlay, it does not render where it should:

In the example of the pictures I've gone from a map size of 1920w 956h to 2012w 1002h, so just a slight expansion of the map size.

I do believe that I know the reason for this, but I don't know how to handle it. I assume the problem is related to the fact that when map changes size - the viewport doesn't, so the getCoordinateFromPixel() is returning "wrong" coordinates for the new map size.

I verified this "assumption" by getting the line/geometry of the measure, and printing the geometry.getLastCoordinate() in pixel coordinates both before and after the map.setSize() method:

var size = [1920, 956];
var coord = getPixelFromCoordinate(geometry.getLastCoordinate())
coord with original map size 1920x956 = [1169,439]

//increase map size and set new resolution
var printSize = [2012, 1002];
map.setSize(printSize);
var scaling = Math.min(width / size[0], height / size[1]);
map.getView().setResolution(viewResolution / scaling);

//print pixel coordinates again
coord = getPixelFromCoordinate(geometry.getLastCoordinate())
coord with new map size 2012x1002 = [1169,439]

So despite the map size was increased, the pixel position is still the same for the feature - which I believe comes from the fact that the viewport did not change.

The API documentation for getPixelFromCoordinate() states "A pixel position in the map viewport." - so this is correct behavior I assume, but what can I do here to "translate" the position correctly according to the new map size?

If there is a simple apparent solution here I will apologize upfront for spending your time to read this, but I have made multiple attempts and not yet been able to find a solution to this.

Kind Regards,

Eivind Rønnevik


texts:
I found a very nice tutorial at: http://sproke.blogspot.com/2012/02/mapping-software-for-84.html, which has some "unix-like scripts", which should be use with QGIS. How can I copy and paste these scripts into MSYS?? Starting MSYS on my win-machine there is no menu-bar and strg+V doesn't work

PS. I search the qgi-manual for vers. 1.73, but found no info. any hints, links ??


texts: St_Transform() - Postgis function Error: column "geom" does not exist

texts:
I am trying to add a new column in the planet_osm_polygon table in order to have a new column with the areas of all the polygons in the database. These are the commands I have used:

psql -U postgres -d gis -c "alter table planet_osm_polygon add column area real"

psql -U postgres -d gis -c "update planet_osm_polygon set area=st_area(st_transform(geom, 2154))

How am I supposed to use ST_Transform() function?


texts: Finding middle point (midpoint) of line in QGIS?

texts: How can I locate the coordinates of the middle point (midpoint) of a line feature in QGIS?


texts: Alternative layperson-friendly names for 'Legend'?

texts:
I am making a mobile mapping app that has a button that displays a legend. Because my legend is initially hidden from view, it is really important the text I put in the button that displays the legend, be understandable to the decent amount of people who dont know what the word "Legend" refers to on its own.

Since it is a mobile app, I only have so much space available, so shorter suggestions (approx. 10 characters) are ideal.

All the alternatives I can think of are equally ambiguous:

-Icon Key
-Marker Key (my data is all points)
-Map key
-Info/ the 'i' info symbol
-the icon seen on many mobile apps, three parallel horizontal bars

Ideas?


texts: How to improve performance of nested search cursors?

texts:
I have three nested SearchCursors. It works correctly but with slow performance. There are some solutions to improve the searchcursor performance like dictionary but i don't know how to use it with conditional statements. The code is as below.How to get rid of nested loops and improve performance?

with arcpy.da.SearchCursor(parameters[1].valueAsText,["Station","Code"],where_clause = sql_exp1) as cur:
    for row in cur:
        if  row[0] != None:
            with s arcpy.da.SearchCursor(parameters[2].valueAsText,["Code","Name","Code1"]) as curt:
                for rowt in curt:
                    if rowt[0] == row[1]:
                        with  arcpy.da.SearchCursor(parameters[0].valueAsText,["Name","Origin","Destination","Code1"]) as curl:
                            for rowl in curl:
                                if rowl[0] == rowt[1] and rowl[3] == rowt[2]:
                                    listnumber.append(rowl[0])
                                    linecodeFinal.append(rowl[3])
            codelinestr = ','.join(str(b) for b in linecodeFinal)


texts:
As newbie with PostGIS, I downloaded Postgresapp for MacOSX. Everything works fine, but I noticed some differences with Stackbuilder which adds to extensions address_normalizer and pgRouting.

Can I add those extensions or are they not necessary?

If necessary how can I add them?


texts: 3d Intersect query compatible with GeoTIFF/DEM?

texts:
When researching this topic I haven't been able to verify if it's possible to make spatial 3D Intersect queries between a vector polygon with z-values and GeoTIFF/DEM file holding the elevation values?

The intersect query I want to run should basically take a planned house (polygon with building height stored as a z-value in the attribute table) and check if it spatially intersects in the z-dimension with height limitations calculated with the plugin advanced viewshed analysis (QGIS).

reading PostGIS documentation on ST_Intersects and ST_3DIntersects rules out the standard ST_intersects query since its only two dimensional. And from what I can understand the ST_3DIntersects query can't work with GeoTIFF's/DEM's? Only TIN-images if SFCGAL backend has been enabled.

  ST_3DIntersects — Returns TRUE if the Geometries "spatially intersect" in 3d - only for points, linestrings, polygons, polyhedral surface (area). With SFCGAL backend enabled also supports TINS

I could probably convert my GeoTIFF's to TIN but I rather not since I have 400 GB of already calculated results from the viewshed analysis.

Is there to your knowledge, any know way to run this kind of intersect query with QGIS/PostGIS?


texts: How to animate raster time series using open-source GIS?

texts: I have a bunch of rasters taken about 10 minutes apart, and I have a nice QGIS project that produces pretty enough maps of them. I'm wondering if there's a way to animate this raster series, with a color map and overlaid with (static) vector data, using open-source GIS. A QGIS solution would be ideal since I already have the project set up there, but other suggestions are welcome as well.


texts: Leaflet Routing Machine: how to initialise map without waypoints loaded?

texts:
I'm starting to take my first steps using Leaflet Routing Machine. Everything works fine if I use load waypoints when I create my L.Routing.control and add it to the map. Now I'd like to start my map without waypoints and after starting to editing points in my map. I've tried to use:

waypoints: [null]

This doesn't work: a blank map appears.

Here's my complete code:

<!DOCTYPE html>

        <meta charset="utf-8"/>
        <title>My Leaflet Routing Machine Example</title>
        <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css"/>
        <link rel="stylesheet" href="LeafletRoutingMachine/dist/leaflet-routing-machine.css"/>
        <link rel="stylesheet" href="http://netdna.bootstrapcdn.com/font-awesome/4.0.3/css/font-awesome.css"/>
        <link rel="stylesheet" href="css/easy-button.css"/>

In this case this is the result:

Otherwise if I use this Routing.control definition:

var control = L.Routing.control({
    //   waypoints: [null],
    waypoints: [
        L.latLng(44.91221, 7.671685),
        L.latLng(44.907852, 7.673789)
    ],
    routeWhileDragging: true,
    show: true,
    language: 'it',
    geocoder: L.Control.Geocoder.nominatim(),
    autoRoute: true
}).addTo(map);

The result is this:

Any suggestions, examples, or jsfiddles to help me?


texts: Count features in polygon by using aggregate (virtual field)

texts:
I have a points layer and a polygon layer. The points are placed randomly. Some are within a polygon, some are outside. I would like to write the amount of points inside each polygon as virtual field into it's attribute table. I know, the aggregate function should do the job, but it don't:

aggregate('dots','count','id',within($geometry,geometry(@parent)))

I get 0 as result in every (multi-)polygon.


texts: Statistical analysis of polygons QGIS : find areas of statistical difference

texts:
I have polygons for the country divided into E.D. areas.

These polygons illustrate the average annual rainfall in their ED area.

Is there any way to statistically tell if one polygon area has a statistically higher or lower volume of rain compared to its surrounding areas?

or if I extract points from these polygons and make a smooth surface raster via interpolation - is there a way I can I assess my layer for statistically higher or lower rainfall ?


texts: How to count K Means iterations on Earth Engine?

texts:
I've used ee.Clusterer.wekaKMeans on a Sentinel 2 image and I was hoping to know if there's a way to see how many iterations were done after the clustering is complete.
I know that one of the parameters on ee.Clusterer.wekaKMeans is 'maxIterations' but I set it to null, so how can I see how many iterations were done afterall?


texts: Changing the URL ID/Handle of an ArcGIS Online Web App

texts:
My organization has created a series of web apps in ArcGIS Online (using Software as a Service). We plan to create a series of regularly scheduled products, and we want the URLs to be more intuitive. For example:

Rather than:

https://organization.maps.arcgis.com/apps/webappviewer/index.html?id=6647cf0dd0a959

We would want:

https://organization.maps.arcgis.com/apps/webappviewer/index.html?id=Product22DEC2016

Notice the updated "id=" value.

I have tried updating the URL in both the general settings in ArcGIS Online (below):

I have also gone into https://ago-assistant.esri.com/ to update the URL.

Both methods just change the reference, and cause (essentially) a 404 error where the product can no longer be found.

Is there a way to actually alter the URL Handle or the ArcGIS Online App ID?


texts: Problems with projections in QGIS 2.2 on a Windows Vista?

texts:
I use QGIS 2.2 on Windows 8, and the same QGIS 2.2 has been installed on a computer with Windows Vista.

On the Windows 8 everything works just fine, but something is wrong on the Windows Vista.
All the projections that I can choose from on the Windows 8 isn't accessable from the Vista computer?

The national projection that all our data are in, can not be chosen on the Vista?
No problem on the Windows 8.

What may we have done wrong? We tried to reinstall QGIS on the Windows Vista, but it didn't help.

Something else we can try?

EDIT:
I tried installing QGIS in a folder without blanks, as suggested in a comment below. But the problem is still there. I still can't use the right projection on the Vista.
The answer in the other question is too complicated for me.

From the postinstall.log (partly in Swedish):

"C:\QGIS_Valmiera\QGIS" „r inte ett internt kommando, externt kommando,
program eller kommandofil.
OSGEO4W home is C:\QGIS_Valmiera\QGIS Valmiera
Synchronizing CRS database with GDAL/PROJ definitions.
Load srs db from: C:/QGIS_Valmiera/QGIS Valmiera/apps/qgis/./resources/srs.db
QFSFileEngine::open: No file name specified
QFSFileEngine::open: No file name specified
QFSFileEngine::open: No file name specified
QFSFileEngine::open: No file name specified
QFSFileEngine::open: No file name specified
Loading epsg.wkt
0 WKTs loaded
CRS update (inserted:0 updated:0 deleted:3631 errors:0)
No CRS updates were necessary.


texts: Speed up rendering background vector data in QGIS

texts:
I have a geopackage of ordnance survey background data (field boundaries, rivers, roads, etc). I won't be doing any analysis on this data it will purely be used as a background / basemap.

My problem is that it's very detailed and quite slow to render in QGIS. Is there a way of speeding up the rendering somehow? I don't know anything about creating vector tiles but is that an option?


texts: Raster in QGIS with Python: blank window

texts:
I'm trying to learn to write QGIS Python Applications. To start with something simple, I thought I'd try to write an application to open a raster layer in QGIS. Here is the code:

import os

from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4 import QtGui

app = QtGui.QApplication([])
QgsApplication.setPrefixPath('/usr', True)
QgsApplication.initQgis()

print "Qgis initialised"

window = QtGui.QMainWindow()
window_frame = QtGui.QFrame(window)
window.setCentralWidget(window_frame)
frame_layout = QtGui.QGridLayout(window_frame)

canvas = QgsMapCanvas()
frame_layout.addWidget(canvas)

filename = 'filepath'
fileinfo = QFileInfo(filename)
path = fileinfo.filePath()
basename = fileinfo.baseName()

if fileinfo.exists():
    print "raster exists"

if QgsRasterLayer.isValidRasterFileName(path, "unknown"):
    print "File is a valid raster"
else:
    print "File is not a valid raster"

layer = QgsRasterLayer(path, basename)
QgsMapLayerRegistry.instance().addMapLayer(layer)
canvas_layer = QgsMapCanvasLayer(layer)
canvas.setLayerSet([canvas_layer])
canvas.zoomToFullExtent()

if layer.isValid() is True:
    print "Layer working"
else:
    print "Layer is not working"

window.show()
print "window shown"

QgsApplication.exitQgis()

app.exec_()
print "app executed"

This is the output when I run it:

Qgis initialised
raster exists
File is a valid raster
Layer working
window shown

and a window does open, but it remains blank. The filepath is correct, if I simply write in the bash console:

qgis filepath

qgis opens it correctly. I've tried changing the canvas color:

canvas.setCanvasColor(Qt.yellow)

and it works; it's just the image that doesn't appear.

I don't understand if it's a code problem or a system problem. My Linux distribution is Fedora 23.


texts:
Background:
I'm using GPS + pressure/temperature/humidity sensor (for elevation, BME280) + 10 DoF IMU (BNO055) to collect data about small geographical region (few hundret-ish square meters) by walking along specific, preselected paths, sidewalks and roads.

I intend to take the walks along the same paths several times to have as much (redundant) data to work with later as possible. My goal is to fuse GPS data with pressure data and IMU data to create three-dimensional map of the region.

GPS data rate is 1 data point per second.
Pressure/temperature/humidity data rate is 1 data point per second.
IMU data rate is ~40 data points per seconds.

At this point all the data is stored in three separate text files (NMEA + CSV + CSV) with time stamps for each data point. The time stamps is what currently keeps these data sets "synced".

Question:
Now that I have this set of three files (with new sets added every time I repeat the measurements), what would be the best method to organize the data? Is there any data format that is intended to combine different readings (with different data rates) into one, coherent stream of data? Is such format or method of organization dictated by work flow?
Which work flow would you recommend for working with GPS + pressure + IMU data sets?

I don't want to discard any data points. I don't want to lose raw data's resolution. I want to keep original data's format.

I'm not asking about data fusion itself, only about organizing collected data.

Please let me know if this is not the best place to ask this question, perhaps I should have posted it to Robotics.


texts: rasterize vector creates raster with all cells containing the same value but different color

texts:
I have a vector layer which i am trying to rasterize. It is smaller than the study area. the end goal is to create a layer with euclidean distance from said rasterized vector layer.
When rasterizing, I use another raster layer covering the whole study area for the output extent. I think this should give me the value 1 for "presence" of my vector and 0 for "absence" in the respective locations. However, I get the value 1 for all cells in my raster layer but the color is differentiated according to the vector i rasterized. The black pixels form the correct shape but the surrounding transparent pixels have the same value. How do I get this raster with the transparent pixels set to 0?
I followed the answer to this question: Calculating distance to points in QGIS and added a new column to my attribute table which i populated with 1 before rasterizing. I used this field as my "field to use for burn-in value". I believe this might be where my mistake stems from but if I leave that field empty i get a completely empty raster layer with values that don't make sense (see layer "Rasterized" in screenshot).
I also tried to set the accuracy in the layer styling panel to actual as had been suggested to me on a different question (reclassify raster based on pixel color) but to no avail.
Thanks in advance!


texts: Is there a standard projection for Sub-Saharan Africa?

texts:
With the help of d3.js, TopoJSON, and data from Glottolog I have made this map of Bantu languages that exist in Sub-Saharan Africa.

I've used the standard Mercator projection so far. I know that there are always tradeoffs between different projection properties but I'd like to know if there is some standard projection that is recommended for Sub-Saharan Africa other than Mercator.


texts: Understanding resolution in vector maps?

texts:
I'm developing a vector maps with detailed features such as road edges, cross walks, traffic lights , street lights etc., by using 3D points (x,y,z) without geo-referencing, generated from stereo camera. But I don't know the resolution of the maps I'm generating.

How I can measure the resolution of my vector maps I'm developing?


texts: Looking for MapGuide 2.1.0 web extensions / web studio - where am I going wrong?

texts:
I am trying to get to grips with MapGuide Open Source 2.1.0, and the first stumbling block is the Windows installer (!).

I got the main binary from the MapGuide download page which has the description 'executable that will install the MapGuide Server and MapGuide Web Extensions'. From what I can gather, previous versions of MapGuide (e.g. 2.0.2) had two separate installers - one for the server and another for the web extensions, but the 2.1.0 installer is supposed to bundle the two.

This being the case, why can't I see a web studio directory inside C:\Program Files\OSGeo\MapGuide\Web\www ? It seems like web studio is not installed.

I tried installing the 2.0.2 web studio from its separate installer but that seems to be causing problems.

Can anyone tell me where I'm going wrong?

Any help much appreciated.


texts: Conversion of MLS/AURA RELATIVE HUMIDITY PRODUCT to correctly referenced GeoTIFF

texts:
I have downloaded the following product

https://cmr.earthdata.nasa.gov/search/concepts/C1251101764-GES_DISC.html

and it gives me a .he5file.
Can anyone help me in correctly geo-referencing the file into a corresponding GeoTiff file, so it correctly projects itself in something like QGIS?

There is also an associated .xml file with it available to download.


texts: If converting .xyz to .las using e.g. laspy, then how to specify offset and scale?

texts:
I have data from Lidar in .xyz format.

I read that it could be beneficial to convert it to .las. So that one can use e.g. laspy -library.

For example (code copied from an example):

import laspy
import numpy as np

hdr = laspy.header.Header()

outfile = laspy.file.File("output.las", mode="w", header=hdr)
allx = np.array([1.000, 2.000, 3.000, 3.000]) # Four Points
ally = np.array([0.000, 0.000, 0.000, 3.000])
allz = np.array([10.000, 10.000, 11.000, 11.000])

xmin = np.floor(np.min(allx))
ymin = np.floor(np.min(ally))
zmin = np.floor(np.min(allz))

outfile.header.offset = [xmin,ymin,zmin] # how to know what this should be?
outfile.header.scale = [0.001,0.001,0.001] # how to know what this should be?

outfile.x = allx
outfile.y = ally
outfile.z = allz

outfile.close()

However, I don't have access to any "header data" or metadata.

How can I know what to specify offset and scale as?


texts:
For using liblwgeom function in spatialite, I need to compile the library. I found that it is included in the PostGIS code.

My question is, how to compile the library without compiling the whole PostGIS?


texts: Iteratively rename feature classes using geodatabase names with Python

texts:
In a folder I have a group of file geodatabases that have unique names from one another. Within each geodatabase is a polyline feature class that is not uniquely named -  they are all called 'Elev_Contour'. I want to use Python to iterate through the geodatabases and append the geodatabase name onto the end of the feature class names to make them unique. After that is complete, I want to copy the feature classes into one new geodatabase.

I've used da.walk to build a list of all the feature classes, but I am struggling with how to access just the geodatabase name when iteratively renaming the feature classes.

import os
import arcpy
from arcpy import env
env.workspace = r"b:\contours"
fclist = []
for dirpath, dirnames, filenames in arcpy.da.Walk(env.workspace,datatype="FeatureClass"):
    for filename in filenames:
        fclist.append(os.path.join(dirpath, filename))

for fc in fclist:
    print fc
b:\contours\Elev_321168_Albany_E_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321169_Albany_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321240_Binghamton_E_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321241_Binghamton_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321259_Boston_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321278_Buffalo_E_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321279_Buffalo_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321420_Elmira_E_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321421_Elmira_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321482_Glens_Falls_E_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321483_Glens_Falls_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321516_Harrisburg_E_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321517_Harrisburg_W_1X1.gdb\Elevation\Elev_Contour
b:\contours\Elev_321600_Kingston_E_1X1.gdb\Elevation\Elev_Contour


texts: Changing symbology from simple marker into a vector (print of a building)

texts:
I've got a question about changing the symbology of a certain rule-based layer. See below as to the current simple marker:

As you can see the house is given this grey color because of the added 'topo' vector layer. The question I have now is, is it possible for me to make the red marker into the vector: I want the whole building to be red. I tried the different symbology settings (raster image, vector field marker etc.) but no luck. It is also (in my understanding) not possible for the vector layer to be rule based because it does not contain the information to be filtered. Maybe combine the vector layer within the other layer that contains the rules?
I can answer the follow-up questions.


texts: Is there a way to manage spatial files in QGIS as in ArcGIS Manager?

texts: I want to copy and paste shapefiles from one folder to another. Is there a way to do it in QGIS as in ARCGIS Manager and not from the Windows Explorer?


texts: Python code in ArcGIS field calculator error?

texts:
I've tried writing a code block using the Python Parser of the ArcGIS Field Calculator. Here's what I am trying to achieve:
There's a field in my attribute table with values for some, but not all of the polygons in the feature class. For those which have a Null Value I need to calculate the values using two other fields and a simple fraction. So the logic should be like this:
If Field1 = Null, new value is Field2 / Field3, else Field1=Field1

I've written a script, but it won't work. It does compute, but there are either no results (meaning that the Null value stays like it is) or, in some cases for some reason, the new value is 1.
This is the script, can anybody see the error?


texts: Splitting worldwide roads into multiple road datasets with less than threshold number of features?

texts:
I want to split a worldwide 'road dataset' into multiple grids/blocks-extents, depending on the amount of features resides in that block-extent. If the amount of features(inside that block) is still above the count of features I want, then it must split the first block into too pieces and so on..., until it has the right amount of features inside that block-extent before I can use an function like PolygonToRaster and after that the mosaic function.

I came up with this route/possible solution, because i have a lot of errors while trying to raster-rise worldwide road datasets.

  ExecuteErrors:ERROR 999998 and 999999: Unexpected Error.

  Failed to execute(PolygonToRaster).

So, if anyone have some experience with this issue, I would gladly know how to solve this issue!


texts: Updating attribute field with random string value

texts:
I have a shape layer in qgis with blank attribute table. I want to fill attribute table with false data randomly. I successfully update the field with integer values using 'rand' expression. I am stuck at 'how to randomly fill string values.'

Suppose I have shape with drain 500 features. I want to enter values in drain_name field with values "ravi", "Sangpo","ganga",... etc.

Which expression should i use to update the drain name with random string value.


texts: How to change Map Scale programmaticaly in C# WinForm Project

texts:
I have created form with axMapControl and comboBox with values of possible scales (e.g. 1:1 000 000, 1: 500 000, etc).

How can I change current map scale programmaticaly by using value from selectedItem taken from Combobox?


texts: How to get feature id for particular feature after feature ismodify in OpenLayers?

texts:
I am using modify control in OpenLayers 2.12. After feature is modified I am running the function below:

building.events.register('afterfeaturemodified',building, onchanges);

Where building is my WFS layer. Now in "onchanges" function runs and now I want my WFS feature id/name, etc.

My function is as follows:

function onchanges(ev)
        {
            var polygonch=ev.feature.geometry;
        }

Now in this function I want feature id/name on which modification is done. Or other method to get feature info on which editing/modification is done.


texts: How to calculate custom fields in ArcGIS Reports?

texts:
I have a custom report in ArcGIS that summarizes my fields using basic descriptive statistics. For example, here is a summary field that calculates the maximum value in the group for a particular field:

What I would like to do is calculate the range (|max-min|). However, this function is not included in the summary statistics. I think that I can use the "Fields>Calculated" menu to add custom fields. Here is a screenshot of that:

Does anyone have an idea of what has to go in as the "Formula" for the field calculation, and how do I fill out the element properties for the field I would like to populate with this calculation:


texts: Assigning points to all bends on polyline in QGIS 3.6

texts:
I have a series of polylines representing water pipelines that I want to clean up.  I want to split the pipeline polyline at key locations where there is a fitting such as a stop valve or hydrant.

Where there is a formed bend I want to add a point and split the pipeline at the bend. Bends come in set angles (say 5, 11.25, 25, 45 and 90 degrees). The pipeline may have smaller bends, say up to 3 degrees, but I want to ignore these.

Is there a way to assign a point at bends on a polyline and assign the angle of the polyline as an attribute? I can select and remove all of the small-angle points.


texts: Loading print layout using an existing template (.qpt) with PyQgis3

texts:
I'm trying to load a print layout (.qpt) which I saved earlier via PyQgis.

I found this earlier entry, but it's for QGIS 2.x and thus doesn't work for QGIS3:
Programmatically load composer from template and generate atlas using pyQgis

Does anyone know how to do this in QGIS3?


texts: How to update PostGIS 2.1 to use `ST_DelaunayTriangles()`

texts:
The error "lwgeom_delaunay_triangulation: GEOS 3.4 or higher required"
arrives on my terminal, running psql with

SELECT PostGIS_Lib_Version();
--   2.1.4

SELECT postgis_version();
-- 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

and  when use use ST_DelaunayTriangles() (ex. in this script).

So, what is the problem, how re-install PostGis or workaround it?

EDIT for @janechii suggestion,

 select PostGIS_full_version();

POSTGIS="2.1.4 r12966" GEOS="3.3.3-CAPI-1.7.4"
PROJ="Rel. 4.7.1, 23 September 2009" GDAL="GDAL 1.9.0, released 2011/12/29"
LIBXML="2.8.0" LIBJSON="UNKNOWN" TOPOLOGY RASTER

PostGIS was installed all by standard apt-get in my Debian Stable this week, and now redo (try reinstall GEOS to use 3.4 not 3.3) with

sudo apt-get -y install postgis postgresql-9.3  libgeos-3.4 libgeos-c1 libgeos-dev

but no effect.


texts: Performing supervised classification with Landsat 7 in Google Earth Engine

texts:
I am trying to perform landcover classification of Shanghai using JavaScript in Google Earth Engine. I am a complete beginner at this and I have been going through many tutorials trying to figure things out. My problem is that many of these tutorials are aimed at Landsat 8 but I need to use Landsat 7 and it covers the time period I need.

I have made my training points and now need to create the training data, so far my code is as follows:

//filter image collection
/**
* Function to mask clouds based on the pixel_qa band of Landsat SR data.
* @param {ee.Image} image Input Landsat SR image
 * @return {ee.Image} Cloudmasked Landsat image
 */
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
};

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filterDate('2013-01-01', '2015-12-31')
                  .map(cloudMaskL457);

var visParams = {
  bands: ['B3', 'B2', 'B1'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};
Map.setCenter(121.51520494843749,31.22304570235914, 10);
Map.addLayer(dataset.median(), visParams);

//merge the feature collections
var newfc = Urban.merge(Water).merge(Vegetation);
//print (newfc);

//create training data by overlaying the training points on the image
var bands = ['B2', 'B3', 'B4', 'B5', 'B6_VCID_1', 'B6_VCID_2', 'B7'];
var training = image.select(bands).sampleRegions({
    collection: newfc,
    properties: ['landcover'],
    scale: 30
});

When I run this it says that image is not defined (after var training), which is true, I have not defined it, but when I try other things (for instance cloudMaskL457.select) it says that cloudMaskL457.select is not a function.

I don't know how to get past this, please help!


texts: Using Python script as precondition in ArcGIS ModelBuilder?

texts:
I have a model that has two preconditions, in which I need to figure out if there is a way to prioritize one precondition to 'run' prior to the second. Is anyone aware of any built-in ability to do so with Model Builder or that can suggest a work around?

For instance, if I try to 'link' one precondition to another, I get the error "Can only connect a variable to a process", this is because my 'output' has a precondition that a python script needs to be run against it first, and then that output has the second precondition executed (but only once the original output has been altered by the Python script).


texts: How to pull lat/lng from a map canvas and update a database

texts:
So if i want to pull out lat/lngs after user clicks on a map canvas (lat/lngs are pulled using any javascript API (leaflet, wax etc.) and add this to a postgresql database using a query, how do i do this?

I came across:

https://github.com/commandprompt/postgres-js

But how do you generally add values to a db on the fly?


texts: Displaying a value as a symbol in Carto?

texts: I am working in Carto and have a geography (NYC Community Districts) and values associated with each community district (# of units built). I is possible to display those values as graduated symbols on a map?


texts: Need to Add/Remove .png Legend with Leaflet Layers Control?

texts:
Supporting Links

Website Location:
http://bit.ly/2fXCiTH
.JS File Location:
http://bit.ly/2g56W0U

Trying to accomplish

I would like to have a specific .png legend associated to each individual layer.

Problems so far

I have not been able to get code to recognize the layers being fired

My non-working code so far

var temperaturesLegend = L.control({
    position: 'bottomright'
});
var xwLegend = L.control({
    position: 'bottomright'
});

temperaturesLegend.onAdd = function(map) {
    var div = L.DomUtil.create('div', 'info legend');
    div.innerHTML +=
        '';
    return div;
};

//temperaturesLegend.addTo(map);

xwLegend.onAdd = function(map) {
    var div = L.DomUtil.create('div', 'info legend');
    div.innerHTML +=
        '';
    return div;
};

// count plus 1
others = 1;

// Add and remove layers
map.on('overlayadd', function(eventLayer) {
    // Switch to the temperaturesLegend...
    if (eventLayer.name === 'temperatures') {
        // xwlegend.removeFrom(map);
        temperaturesLegend.addTo(map);
    } else { // Or switch to the xwlegend...
        if (others < 1) {
            xwLegend.addTo(map);
        }
        others++; // temperaturesLegend removeFrom(map);
        //xwlegend.addTo(map);
    }
});

// Add and remove layers
map.on('overlayremove', function(eventLayer) {
    temperatures = eventLayer;
    //console.log(eventLayer.name);
    // Switch to the xwlegend...
    if (eventLayer.name === 'temperatures') {
        this.removeControl(temperaturesLegend);

    } else { // Or switch to the xwlegend...
        others--;
        if (others < 1) {
            this.removeControl(xwLegend);
        }

    }
});


texts:
I have the following traceback and I would like some help to understand/solve it.
Can anyone explain to me what to look here?

6: stop("cannot use this formula, probably because it is not vectorized")
5: .overlayList(x, fun = fun, filename = filename, forcefun = forcefun,
       ...)
4: .local(x, y, ...)
3: raster::overlay(x, fun = tslp, ...)
2: raster::overlay(x, fun = tslp, ...)
1: raster.kendall(raster_stack, tau = TRUE, intercept = TRUE, p.value = TRUE,
       z.value = TRUE, confidence = TRUE)

example data:

dates<-as.data.frame(
  as.Date(c('2010-10-1','2009-5-26','2008-5-10','2007-3-13','2006-2-2')))

rst = raster(matrix(1:15,3,4))
raster_stack = stack(rst,rst,rst,rst,rst)
raster_stack <- setZ(raster_stack, dates[,1], "my_dates")


texts: Problem in exporting Annotation data to shp

texts:
when transferring my .dgn file to .shp in Arcgis, I can't export the annotations to .shp !!! I do not know why the outil "export data" is disabled !!!

However, when importing the file .dgn in Arcgis, I can see the annotations but the problem is it can't be exported to .shp

How would I fix this?


texts: Displaying hundreds of thousands points on web map?

texts:
I'm dealing with around 150,000 points on a map. My customers need to see all the points, wthout clustering. I tried Openlayers and Mapbox GL Js, both have the same issue: its memory costs are huge!

I stabilized on Mapbox GL JS 0.15 (max. 1GB RAM during loading, min 700MB RAM) but since it's not working on Chrome (some tiles are freezing) I have to go with the new version, which is nice but has a memory map (for my case) of 1.6GB!

My biggest customer has a computer with 4GB RAM, this is clearly not working...

I know this is my fault to display all that points on the client side, but I don't have any experience of finding a good & free solution, client side or server side.

Example of my map on Mapbox GL JS 0.26 (loading is long as I load my data from dropbox)


texts: Create polylines from point array Python

texts:
I am working on building polylines from an array of points I have created from a csv file and am completely lost in building the lines. My code works as follows: I have the point array working, I need help in creating the polylines from the point array. Please note I am also wanting to add the name into the new polyline feature.

# Reads a csv file to create points of locations
import arcpy
print "==============RESTART===================="

# Set up variables
spreadsheet = "C:\\Users\\Desktop\\RhinoObservations2.csv"    # input file with coordinates
ouputFC = Rhino_Path          # defines the name of the ouput feature class

# Open the CSV file and read the header line. "r" is read the file
observationsFile = open (spreadsheet, "r")
print observationsFile

headerLine = observationsFile.readline()  # reads all lines in CSV??
print headerLine

# Determine the indicies of the columns "X","Y","Rhino"
fieldList = headerLine.split(",")
print fieldList

xIndex = fieldList.index("X")
print xIndex

yIndex = fieldList.index("Y")
print yIndex

rhinoIndex = fieldList.index("Rhino")
print rhinoIndex

# Set up new dictionary where the Rhino information is stored
rhinoDictionary = {}

# Loop through the remaining line of the input file; for each line....
for line in observationsFile.readlines():
   print line

    # Determine the rhino name, X,Y coordintes from that line
    segmentedLine = line.split(",")

    rhino = segmentedLine[rhinoIndex]

    x = segmentedLine[xIndex]

    y = segmentedLine[yIndex]

    # Create an object of the the class Point (defined in the arcpy module) for the X Y coordinate
    coords = arcpy.Point(x,y)

    # If this rhino is not already in the dictionary
    if not rhino in rhinoDictionary:

        # Create a new object of class Array, add the object point, and then put array into dictionary
        # under the name of the Rhino
        coordArray = arcpy.Array()  #setting up an empty list to store points
        coordArray.add(coords)      #list that contains one point
        rhinoDictionary[rhino] = coordArray

    # Else
    else:

        # Get the array object from the dictionary for this Rhino and add the point to the array.
        coordArray = rhinoDictionary[rhino]
        coordArray.add(coords)

# Print out the information in the dictionary
for key in rhinoDictionary:
    print key
    print "This is the point array: " + str(rhinoDictionary[key])
    for p in rhinoDictionary[key]:
        print "(" + str(p.X) + ", " + str(p.Y) + ")"

    # Take the array and create polylines from it.

    print "========"
# Create polylines and add them to the feature class


texts: Exporting to KML file with specifying some symbology (icon styles)

texts:
I created a plugin which allows to export layers from QGIS to KML files, and I want to specify the type of icon of the points in Google Earth.
Some one can tell me if there is a way to do that with QgsVectorFileWriter.writeAsVectorFormat?
For example to make this icone for my point in my layer.
This is my code
layer = QgsProject.instance().mapLayersByName('bpe')[0]
output_layer = r"C:\Users\nassim.mokhtari\Desktop\X.0140253.0.xx_FREE_TCA-SITE-CAB14\03- Livrables\03- Shapes\bpe_test_iface.kml"
name = "Name"
datasource_options = []
datasource_options += ["DescriptionField= "+ name]
datasource_options += ["NameField=" + name]
QgsVectorFileWriter.writeAsVectorFormat(layerRoot, output_layer, "utf-8",layerRoot.crs(), "KML",datasourceOptions=datasource_options)


texts:
Hypothetical situation: I've got a database with 100,000 customer addresses, which are spread out across the continental US. I've also got a year end totals for number of purchases and value of purchases.

I'd like to take that table and create districts based on sales data. To start with, I'd like as many districts as necessary with roughly 1,000 purchases in each district. I'd like to avoid gerrymandering as much as possible. I'd like to do the same thing per $1,000,000 in sales. I'd also like to be able to generate districts quickly and repeatedly with slightly different parameters (2,000 purchases, $500,000 sales, etc).

I anticipate that areas with higher population density would be composed of smaller districts. I'd use each of these iterations to manually define my own districts based on available personnel, etc.

Here's my problem: I know this can be done, but I'm not even sure which search terms I should be using to find some guidance. I've read up a little on election redistricting, but I am having a hard time figuring out how to make that approach work for me. It's not just where the customer/voter is located, but there are also financial components (such as sales volume), and I don't see those kind of considerations used in creating election districts.

What do you call what I'm trying to do? What tools/procedures/processes should I begin to dig into? I'm prefer a solution using QGIS, but I have the entire OSGeo suite. I do not currently have access to any ArcGIS products.


texts: Getting 2147483647 and -2147483647 as min and max output for point to raster, polygon to raster, and feature to raster tools

texts:
No matter which feature to raster tool I use (feature to raster, point to raster, or polygon to raster), I am getting a raster output whose min and max values are 2147483647 and -2147483647, respectively.

Researching the issue came up with nothing much except a Wikipedia article mentioning that this is "the maximum value for a 32-bit signed integer in computing." I was relieved to discover that the number does, in fact, have some rational meaning, albeit one that doesn't make much sense to me...

I am attempting to convert evenly spaced (500 ft) points (and/or 500 ft vector grid) of elevation values to a raster. I'm not interested in interpolation. I simply want the raster grid to reflect the vector grid.

There are 6718 records, whose values, in increments of 100, range from 800 to 1600. The field type is a short integer with a precision of 4. I want to give an output cell size of 500 feet.


texts:
I have multiple points and I need to create polygons that consist specific amount of points.

How could i do this?


texts: Replacing last word in string using ArcGIS field calculator?

texts:
In this table I need to change the last word from field !name! using the text from field !abbreviation!

I was trying to use Python parser to get the text from !name! field without last word and then add !Abbreviation! field. But my expression returns an error:

!Abbreviation! = !name!.split(' ')[:-1]

Can anybody help me?


texts:
Having some issues with the tabulate area spatial analyst tool. Not sure why it's not working, but the output table is missing zones from the input feature class that defines them.

I have a simple raster with 3 values. I'd like to simply find out the amount of cells that lie in each zone of the polygon. The tabulate tool should do this no problem.

But, after the tool runs, the output table only has about 1800 records. The original input zone polygon has around 5000. Any ideas on why this tool is failing so badly??


texts: Use Anaconda's Python distribution within OSGeo4W

texts:
I want to use a Python script that executes Grass commands via os.system(). This of course, has to be executed from a running Grass session. However, when I use the OSGeo4W Grass, it does not find all imported packages as it falls back to its own/default Python distribution. How can I enforce OSGeo4W to use the Anaconda Python distribution?

I found a solution on the internet that makes the OSGeo4W distribution available systemwide, but not the other way around. How is it possible to do that?


texts: PyQgis, Add field to attribute table and Output error

texts:
I'm trying to code a script who calculate id of vineyard lines and vines in these lines.

I code with temporary layers, just my final layer will be save.
But I'm getting stuck on adding new fields, at line 39, I have an OUTPUT error, and QGIS log says no more.
Can you help me ?

This is my code :

##point_layer=vector
##Intervine_distance=number 1.5
##output_layer=output vector

import os, processing
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *

# Project Layer in L93

l93 = processing.runalg("qgis:reprojectlayer", point_layer, "EPSG:2154", None)

# Buffer

Dbuf=(Intervine_distance/2)+0.1

buffer = processing.runalg("qgis:fixeddistancebuffer", l93['OUTPUT'],
Dbuf, 100, False, None )

# Dissolve

dissolve = processing.runalg("qgis:dissolve", buffer['OUTPUT'], True, '',  None)

# Split

split = processing.runalg('qgis:multiparttosingleparts', dissolve['OUTPUT'], None)

#New Fields

rang = processing.runalg('qgis:addfieldtoattributestable', split['OUTPUT'], "Rang", 0, 3, 0, None)

pied = processing.runalg('qgis:addfieldtoattributestable', rang['OUTPUT'], "Pied", 0, 3, 0, None)

# Some steps remaind to code here

# Reproject Layer in WGS84

processing.runalg("qgis:reprojectlayer", pied, "EPSG:4326", output_layer)


texts:
I have an attribute called catdet_sector which is of type text.
The problem is that the values of this attribute appear in parentheses as shown in the attached figure.
I present a personalized form which presents the values, but when the attribute is in parentheses it does not present them.

I tried the replace function but it did not work, the problem is because of the type of object I use in the form. When I use the ratio of values or value map in edit mode, the parentheses appear. On the other hand, when I use text editing the parentheses disappear and the data is displayed correctly


texts: Displaying band names in netcdf file in ArcGIS 10

texts:
I have several netcdf files that contain multiple bands with useful names. When I import into R, the names are intact. However, when I import the netcdf files into ArcGIS 10.2 using the makeNETCDF raster layer tool, the names do not appear and I am left with "band1, band2" etc.

Does anyone know how to assign or import band names for netcdf files?


texts: WFS in QGIS : different boundingbox CRS than Feature CRS

texts:
I want to add a WFS 2.0 layer to my QGIS 3.16 and I have a bounding box issue: QGIS sends a bounding box in the same CRS than the layer CRS. But the server expects the bounding box to be in WGS84, whatever the Feature default SRS.
The GetCapabilities states:
<wfs:FeatureType>
<wfs:Name>LUMI_LUMITHD:Base_prise_BAL</wfs:Name>
<wfs:Title>Base_prise_BAL</wfs:Title>
<wfs:DefaultCRS>urn:ogc:def:crs:EPSG::102582</wfs:DefaultCRS>
<ows:WGS84BoundingBox>
<ows:LowerCorner>-29.22026123106372 39.04255677649548</ows:LowerCorner>
<ows:UpperCorner>70.12090695084257 84.78250333201653</ows:UpperCorner>
</ows:WGS84BoundingBox>
</wfs:FeatureType>

So it seems this WGS84BoundingBox is not implemented in QGIS WFS client. It seems it already is implemented for WMS
Can anybody confirm this or point me to a solution ?
How can I make QGIS report a "standard" wgs84 bounding box ?
Note : I can work around the problem by removing the "Only request features overlapping the view extent" flag, but this leads to performance problems.
Note 2 : When I send a request through the browser with the correct WGS84 BBox, the server will answer with the features so it seams to be a problem on QGIS side.


texts: How to calculate the distance between points (pointlayer)

texts:
I am totally new to QGIS and I have a List of Trackpoints which I loaded in QGIS as a CSV (with the GEOCSV plugin to be able to edit the layer).So I have a Layer with ~1000 Points now.

I like to add a column "distance" in the table of attributes with the distance from the previous point. The Distance should be calculated using latitide and longitude.


texts:
Is it possible to use leaflet-draw-toolbar with angular 2 and ngx-leaflet-draw?

I do not see, and I am not versed in Angular 2, enough, to incorporate the JS library for leaflet-draw-toolbar into my code.


texts:
I'm using QGIS LTR 2.14.7 running on debian.
2 versions of GDAL are installed on my computer (1.10 and 2.0).

1.10 is out of box GDAL version (doesn't support ECW).
2.0 is compiled (by me) to support ECW.

In QGIS I change system environment libraries like it:
GDAL_DATA=/usr/local/share/gdal (I've installed GDAL 2 in this path).

But QGIS continues to use GDAL 1.10 version and not GDAL 2.
Have you an idea to force QGIS to use GDAL 2?

I can't uninstall GDAL 1 because QGIS needs this version to start (version of official repositories).

I've tried "dirty" solution by replacing existing symbolic link:

/usr/lib/libgdal.so.1 -> libgdal.so.1.17.1

by:

/usr/lib/libgdal.so.1 -> /usr/local/lib/libgdal.so.20

This method permit to read ECW files.

But all modifications on a integer field make troubles.
I've read that GDAL 1.0 use 32bits integer and GDAL 2.0 uses 64bits integer.
So QGIS will not update field (type error: type not match with field type).


texts: Making Ridge Line and Valley Line from Contours or DEM using SAGA GIS

texts:
How do I make ridge lines and valley lines using open source GIS Software?
I am totally stuck. I found solution in SAGA GIS - also available in QGIS but it freezes when I run it.
http://www.saga-gis.org/saga_tool_doc/2.3.0/ta_morphometry_24.html
I used default settings. I used this method in separate SAGA app and also in QGIS both same issue.
I have DEM 0,25x0,25 in tiff and also countour lines in shp with 1m step.


texts: Calculation proportion of pixels from 2 raster in qgis

texts:
I am new to Qgis and also R, so I am struggling a bit...
I have two raster layers of same resolution that perfectly align. One raster is DEM, and the other one is a binary forest map (forest = 1 , no forest = NA).
All I have to do is to calculate the proportion of forest per elevation. Like, (how many pixels of forest fall in elevations between 0-50) divided by the (total number of pixels in the DEM that have the value between 0-50). How can I calculate this proportion (forest per elevation) in a automatic way?


texts: Why $area gave me a wrong value QGIS 2.99?

texts:
I want to calculate the area with $area function in QGIS 2.99 but the operation gave me a wrong value.

Here is a picture comparing the results between QGIS 2.99 and QGIS  2.18. The right side is QGIS 2.18 and its OK, I compare the result with ArcGIS.

Here are the projection properties from each one.

QGIS 2.99

QGIS 2.18


texts: Adding transparency to grayscale image using GDAL?

texts:
I have some binary mask image with 255 and 0 values and I want 0 values become transparent.

How can I do this?

Seems it work, but I'm not sure if it the best method:

gdal_translate -a_nodata 0 -of PNG src.png dst.png


texts: How to separate or cut the elevation from bathymetry raster TIFF file with vector layer coastline in QGIS 3.4.3?

texts:
I have a raster GeoTIFF file downloaded from GEBCO 2014 which contains bathymetry and elevation.

I need to split in 2 layers the raster GeoTIFF file (layer: Bathymetry and layer: Elevation) using coastline vector layer.

Can QGIS 3.4.4 do this?


texts: Creating a constant surface at designated height above river

texts:
I am trying to create a polygon (and eventually a raster) surface of river deposits that have since been eroded away. We know that at maximum, deposits got up to 600 ft above the present river. Given that a river slopes as it goes down hill, the answer is not as easy as a raster calculation, or determining X ft above sea level + 600 ft. How does one make a sloping surface as constant elevation above certain points in the river to determine the previous geographic extent of these deposits?

Using ArcGIS 10.5.


texts: QGIS plugin problem, can't reach objects in a Qt-based Stacked Widget

texts:
I am creating a plugin for QGIS but am having trouble with the Qt components. I have used Plugin Builder and Qt designer to create all necessary files for the plugin.
I want to use a QStackedWidget inside my QDialog and I can open the dialog just fine when I don't try to access the "pages" (or widgets) of the stacked widget. But I am getting an error when trying to reach the pages/objects located inside the pages of a stacked widget from my code.

Simple code example:

def __init__(self, iface):

    self.dlg = Dialog()
    self.stackedWidget = self.dlg.stackedWidget
    self.dlg.stackedWidget.page.pushButton.clicked.connect(self.func)

def func(self):
    :code that does some stuff:

Returns the error: AttributeError: 'QStackedWidget' object has no attribute 'page'

Similarly if I try to work with just a stacked widget not inside a QDialog (again opens just fine without the code that causes errors) I get the same kind of error.

Code Example:

def __init__(self, iface):

    self.stackedWidget = StackedWidget()
    self.page = self.stackedWidget.page
    self.page.pushButton.clicked.connect(self.func)

def func(self):
    :code that does some stuff:

Returns the error: AttributeError: 'QWidget' object has no attribute 'pushButton'

I am not sure what I am doing wrong or what is causing the error but will of course provide whatever code/information that might be needed to help me solve this problem.

I am thankful for whatever help you can provide!


texts: Preventing pyramid build when loading jp2000 image file into ArcMap or ArcCatalog?

texts:
Every time I want to bring a jp2000 imagery to ArcMap, it has to run its pyramids before it can display it.

Is there a way around for me to prevent from doing it repeat ?


texts: Apply .median reducer to image collection over roi EXCLUDING images with high cloud cover. GEE

texts:
I'm trying to train an algorithm to detect mangrove forests in several regions around Asia for each year using Landsat 8. When I apply the .median reducer on some image collections over a 1-year timespan, the composite is totally clear and cloud free (see 1st image - code link: https://code.earthengine.google.com/486067b4ca23929c45313eadd0e0cacc)

However, for some regions or years with high cloud cover throughout the year, the .median reduction gives me a composite that is noisy (see 2nd image - code link: https://code.earthengine.google.com/286e97966b48d2ae426b832f024097ee).

Is there a way I can apply the .median reducer only on images with comparatively low cloud cover? Or can anyone here think of another way of getting an accurate, cloud-free, year-long composite for an image collection to do a supervised classification?

Code:

    // filtering the image collection
var image = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .select(['B[1-7]'])
  .filterBounds(roi)
  .filterDate('2015-01-01', '2015-12-31')
  .median());
Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 0.3}, 'landsat');

// // merge feature collections
// var newfc = mangrove.merge(other_veg).merge(water).merge(urban);
// print(newfc);

// // create training data
// var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'];
// var training = image.select(bands).sampleRegions({
//   collection: newfc,
//   properties: ['landcover'],
//   scale: 30
//   });

// print(training);

// //Train the classifier
// var classifier = ee.Classifier.cart().train({
//   features: training,
//   classProperty: 'landcover',
//   inputProperties: bands
// });

// //Run the classification
// var classified = image.select(bands).classify(classifier);

// //Display classification
// Map.centerObject(newfc, 11);
// Map.addLayer(image,
// {bands: ['B4', 'B3', 'B2'], max: 0.3},
// 'Landsat image');
// Map.addLayer(classified,
// {min: 0, max: 3, palette: ['008b04', 'ffc82d', 'ff0000', '1400c2']},
// 'classification');
// Map.addLayer(newfc, {}, 'data_points');

https://code.earthengine.google.com/e03d657c9b91c43056f51a175f0a0ff2


texts: How to export maps with PNG images with transparent background to PDF?

texts:
On the Data View in ArcGIS 10, I have some polygons and one PNG image (wind rose) with transparent background.

When I export what is showing on the screen to PNG, the compass appears in the output. When exporting the same to PDF, it doesn't.

Is there a way to get everything exported to a PDF?


texts: How to determine if a tiff is georeferenced or not?

texts: I received a tiff file from someone, but I'm not sure if it's georeferenced. Is there an easy way to check this?


texts: Create single KML from multiple encoded polylines

texts: I have a number of encoded polylines stored in separate database records.  I can successfully retrieve each one, decode it and create a single KML file.  However, I need to create a KML file that includes all polylines.  What I've discovered is that the syntax for both geojson (my intermediate step) and KML are drastically different when you have more than one line segment.  None of the libraries I've found create MultiLineStrings (http://geojsonlint.com/#sample-MultiLineString) and the tokml library (https://github.com/mapbox/tokml) doesn't seem to support them as an input anyway.


texts: How do you symbolize a feature class layer in ArcMap 10 based on one attribute's value and then another attribute's value?

texts:
I currently have a point feature class that has a couple of attributes that I want to use to symbolize its layer.  One of the attributes is called ExceedsTolerance, which either has "Y" or "N" for its value.  The other attribute is called Measure, which has the values -2, -1, or anything from 0 and greater.

I know that I can create a layer definition that has different symbols for unique values for each attribute.  For instance, I can have a black box for Measure when the value is -2, an orange circle for when the value is -1, and a blue circle for when it is greater than or equal to 0.  The same idea goes for ExceedsTolerance.

What I would like, however, is to have it where it shows a certain symbol when ExceedsTolerance is "Y" no matter what the Measure is, but when ExceedsTolerance is "N" to show whatever I have set for Measure.  I can somewhat do this with the Multiple Attributes section of the Symbology tab in the Layer Properties, but I can't do it with the case where Measure is greater than 0 and ExceedsTolerance is "N".  How would this work in ArcMap and is it possible?


texts:
I have a timer script which shows the time of system in Description tag

<script type='text/javascript'>
setInterval('tick()',50);   //tick is a plugin and 50 is in ms delay time
 function tick()
 {
   document.getElementById('clock').innerHTML=new Date();
 }
</script>

I want that if one placemark is placed in India , then same script should show the current time of India and if one placemark is in America , it should show America Time only regardless of machine in which KML has been loaded.

How can I do this (how can I sync all times)?


texts: Exporting ArcSDE feature classes into file geodatabase with ArcPy?

texts:
I would like to use Windows task scheduler to automatically (monthly) create a new file geodatabase and then export feature classes from SDE into the new file geodatabase.  Im using this code, however I get an exucute error that states that the geodatabase already exists.

# Name: CreateGDB.py
# Description: Create a file GDB
# Import system modules
import arcpy
import datetime

out_folder_path = "w:/GISDATA/Database_Backups"
nowstart = datetime.datetime.now()
YearMonthDay = nowstart.strftime("%Y_%m_%d")
out_name = "TEST_WSdata_" + YearMonthDay + ".gdb"

# Execute CreateFileGDB
arcpy.CreateFileGDB_management(out_folder_path, out_name)

arcpy.env.workspace = 'C:\Users\cbgibson\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\TEST_WSdata.sde'
# Local variables:
waterFC = ['WATERLINES', 'GRAVITY']

for fc in waterFC:
# Process: Feature Class to Geodatabase (multiple)
    arcpy.FeatureClassToGeodatabase_conversion(fc, out_name)


texts: Using spatial points to extract values from large raster in R

texts:
I am working with location data for 22 different individuals (id). The data contains an id column, and coordinates (UTM_x and UTM_y). For all individuals combined, there is a total of 991,099 different locations. I am trying to extract raster values (1x1m vegetation classification) for each point location and I am having issues with the speed of the extraction (it takes an extremely long time) and memory issues.
Here is what I have done so far;

First I create a spatial points dataframe from the UTM coords and project them to the correct CRS;

  coordinates(data) <- c("UTM_lon", "UTM_lat")
  data@proj4string <- CRS("+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

Then I load the raster file

  veg.r <- raster("C:/Users/veg.ras.tif")

The raster file was projected in Arc. Checking to make sure the projections are the same;

  proj4string(data) == proj4string(veg.r)
  1 TRUE

Here are the details of the raster;

  veg.r
  class       : RasterLayer
  dimensions  : 81299, 87251, 7093419049  (nrow, ncol, ncell)
  resolution  : 1, 1  (x, y)
  extent      : 606777.517, 694028.517, 4751626.24, 4832925.24  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
  data source : C:\Users\veg.ras.tif
  names       : veg.ras
  values      : 1, 11  (min, max)

Now extract the raster cell values for each point;

  ext <- extract(veg.r, data, df=TRUE)

I have waited 24+ hours for the extraction with no results. I know there isn't an issue with the actual code, because I can
perform this function with smaller subsets of the data.

I have tried using a multicore approach, as suggested HERE, with the code below;

  library(snowfall)
  data.sp <- SpatialPoints(data, proj4string = CRS("+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))

Now, create a R cluster using all the machine cores minus one

  sfInit(parallel=TRUE, cpus=parallel:::detectCores()-1)
  sfLibrary(raster)
  Library raster loaded.
  Library raster loaded in cluster.

  sfLibrary(sp)
  Library sp loaded.
  Library sp loaded in cluster.

  data.df <- sfSapply(veg.r, extract, y=data.sp)
  Error: cannot allocate vector of size 26.4 Gb
  sfStop()

Stopping cluster

As you can see, I get an error due to memory issues.

Are there any suggestions on why the "multicore approach" is not working?


texts: Import SketchUp created KML contour lines model to QGIS

texts:
There's a neat trick to create contour lines based on the GE 3D terrain using Sketch-Up.
The resulting KML however is a model

I have looked at ways of importing KML layers to QGIS, like it was suggested here, yet it doesn't accept this kind of kml file.

So, is there a way to import such a thing into QGIS (or ArcMap for that matter)? and if not can such a model be split into its separate contour lines while retaining the z elevation attribute?

Here's the model itself:


texts: Merging rasters with rasterio in blocks to avoid MemoryError?

texts:
I am trying to merge two rasters with rasterio.

However, when I run dest_array, out_transform = merge(sources), I receive the following error:

MemoryError: Unable to allocate array with shape (1, 56160, 85196) and data type float64

I know rasterio can read and write files in blocks.

Unfortunately, I cannot seem to find a way to do it when calling merge.

Is there a way to use blocks with merge or another way I can merge rasters without using the merge method?

I know how to read and write data in blocks, but I have no idea how to do that when using the merge method.

Here is the code I am using until the line that calls merge and that gives me the error:

import rasterio
from rasterio.merge import merge
import glob
import os

def merge_raster_in_folder(src_dir, dest_dir):
    """Merge all the raster (TIFF extension) within the given directory"""

    # create list of raster files in source folder
    rasters = []
    os.chdir(src_dir)
    for file in glob.glob("*.tif"):
        rasters.append(file)

    # create list of raster objects from the list of raster filename
    sources = [rasterio.open(raster) for raster in rasters]

    # create array representing all source rasters mosaicked together
    dest_array, out_transform = merge(sources)


texts:
I want to erase the memory register that stores the variables that are created while running python code (in QGIS).  I tried using qgs.exitQgis() based on pyqgis manual. But I get an error saying NameError: name 'Qgs' is not defined.

Edit: Realized that qgs.exitQgis() is not usable for my purpose since my aim is not to have a stand alone pycode that uses QGIS library, moreover I noticed that running the command qgs = QgsApplication([], False), in order to define the variable qgs is making my QGIS2.14 interface unresponsive!

Apparently, I need to have some QGIS independent way to clear memory register.


texts: I'm using ArcMap 10.2.2 in a remote server and I have added the online "Imagery with Labels" base map to my map (first picture). When I turn it off the labels (and the light dashed line) remain visible though, and even by removing the basemap from the table of contents, they don't disappear (second picture). Why?


texts:
Platform:
Windows 10
QGis 2.18.14
QT Desiner 4.8.5
PyQT4
I'm creating a custom form using QT Designer to be used in a layer along with a routine in Python.
Just for testing I added a QgsColorButton element on the form.
The idea would be to select the color in the color options box of this button and the selected color would be saved in the attribute table of the layer's feature.
I tried to do a little routine in Python only to check if I could access the properties of the QgsColorButton when the form was opened, but can not even get through this phase.
Any tips to solve this problem?

In the console the code seems to be simple, but I have difficulty understanding with accessing the QgsColorButton widget of the form created in QT Designer.
This would be the form template with multiple QgsColorButton.

I need to access the QgsColorButton with something like the code below:
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *

def form_open (Dialog, layer_id, feature_id):
     global form
     form = Dialog
     btn = form.findChild (QgsColorButton, "mColorButton")
     btn.setColor (QColor (255, 0, 0, 127))
     layer = iface.activeLayer ()
     provider = layer.dataProvider ()
     feature = layer.selectedFeatures ()
     feature ['Color_01'] = btn.color ()
     layer.startEditing ()
     layer.updateFeature (feature)
     layer.commitChanges ()

The code above seems to be simple, but does not work as you like.
A second part code would be the QgsColorButton load the color saved in the attribute table when the form was opened.
I have no idea how to make this work!

@Francisco Raga, perfect! The only change I need to make is to change the LoadColor method.
def LoadColor ():
    color_value = ToolFeature ["Color"]
    if not isinstance (color_value, QPyNullVariant):
    mColorButton.setColor(color_value)
    return

But this returned the following error.
Traceback (most recent call last):
   File "", line 1, in
   File "", line 22, in UpdateColorForm
   File "", line 29, in LoadColor
AttributeError: 'QPushButton' object has no attribute 'setColor'

When I enter the command below in the Python console of the QGIS
dir (QgsColorButton) returns all methods of QgsColorButton and the setColor method appears in the list.
I compiled the ui template to check which methods and arguments I want to return to the method
self.mColorButton.setColor (QtGui.QColor (0, 0, 255))

I do not know what can be done, because in the QGIS console the setColor method is recognized and it runs the code perfectly.
I tried to use the code below, but nothing happens. No error message is displayed and the background of the mColorButton button does not change color with the value loaded from the "Color" attribute of the feature.
def LoadColor():
    color_value = ToolFeature["Color"]
    if not isinstance(color_value, QPyNullVariant):
        style = "QWidget { background-color: %s; }" % color_value
        mColorButton.setAutoFillBackground(True)
        mColorButton.setStyleSheet(style)
        mColorButton.update()
    return


texts:
I have a GeoPandas DataFrame (loaded from a shapefile) containing a number of polygons which looks like this:

                                            geometry  raster_val
0  POLYGON ((8.56799 47.80800, 8.56799 47.77200, ...         1.0
1  POLYGON ((6.73199 47.12400, 6.73199 47.08800, ...         1.0
2  POLYGON ((10.36799 47.01600, 10.36799 46.98000...         1.0
3  POLYGON ((8.60399 46.51200, 8.60399 46.47600, ...         1.0

Is there any elegant way to print/plot those polygons into a .jpg file?


texts:
I am using the following example, to select a vector layer feature on click. The example demonstrates Highlighting and selecting the layer. I was able to successfully implement the Highlighting of the layer but I am not able to select the layer.

This is my code

var highlightCtrl = new OpenLayers.Control.SelectFeature([vector_layer,QuebecSpecies], {
                hover: true,
                highlightOnly: true,
                renderIntent: "temporary",
                eventListeners: {
                    beforefeaturehighlighted: report,
                    featurehighlighted: report,
                    featureunhighlighted: report
                }
            });

            var selectCtrl = new OpenLayers.Control.SelectFeature([vector_layer],
                {clickout: true}
            );

map.addControl(highlightCtrl);
map.addControl(selectCtrl);

highlightCtrl.activate();
selectCtrl.activate();

Where am I falling short?


texts: Polygonize.py: AttributeError: 'NoneType' object has no attribute 'CreateDataSource'

texts:
Given a topographic raster tif of an european country where elevations are simplyfied to either 0, either 200.

How to create a polygon ?

Comment: With the following polygonize.py command:

gdal_polygonize.py level200.tif -f "ESRI_Shapefile" level200.shp

I get to following error:

Creating output level200.shp of format ESRI_Shapefile.
Traceback (most recent call last):
  File "/usr/bin/gdal_polygonize.py", line 164, in <module>
    dst_ds = drv.CreateDataSource( dst_filename )
AttributeError: 'NoneType' object has no attribute 'CreateDataSource'


texts: pyqgis traversing subgroups and group layers

texts:
I am trying to find the name of subgroups using pyqgis. I have been trying to use Playing with subgroups

I am however unable to traverse and get the names of the subgroup. I have tried the code on this article.

for child in root.children():
  if isinstance(child, QgsLayerTreeGroup):
    print "- group: " + child.name()
    if child.name() == "test": #to check subgroups within test group
        for child in root.children(): #which I think i am making a mistake here because the root is supposed to contain all nodes


texts: How to update the items of qCombox in custom form?

texts:
For my project I use a custom form with multiple value relation comboboxes.
This works quite fine and I can access the widgets via my script.

Init the two comboboxes which are correctly filled

cb_type = mydialog.findChild(QComboBox, "type")

cb_subtype = mydialog.findChild(QComboBox, "subtype")

Watch first one to change

cb_type.currentIndexChanged.connect(changeType)

function called when combobox1 is changed

def changeType():
    change filter of value relation of combobox 2 - works
    mylayer.valueRelation(2).mFilterExpression = "\"somefielname\"=\'some attribute\'"

HOW TO REFRESH THE VALUES HERE?

Now the problem is that the attributes are not updated until I close and reopen the form.
Mybe someone can help me to solve the problem, that the list is refreshed automatically!


texts: Cross-browser marker animation in OpenLayers

texts:
Currently I'm trying to implement animated markers in OpenLayers.  I'm using PNGs for the markers, set by the externalGraphic style property on an OpenLayers.Layer.Vector.  I need to be able to make the graphic pulse, grow, blink, etc. which means controlling its size and position.

With the SVG renderer I can get a reference to the image element by using the Geometry's id property.  However with the VML renderer (IE8 doesn't support SVG) I can't seem to find a way to get a reference to the rect element, which is what (I think) I need to manipulate to do the animation.

Am I thinking about this wrong?  Can I change the feature's style on the fly?  Should I just use animated GIFs (ugh)?  OpenLayers.Layer.Markers?  Something else?

What if I could get away with not supporting IE8 -- is working with the image a good way to go?

Edit: I also need to support selection and clustering of the markers...


texts: How to obtain elevation points for a road using contour?

texts:
I have contour lines, & roads centerlines feature classes, I need to get points with elevation in the starts, intersections, & the ends of the roads, how can I get this information? I used the "Intersect" approximity analysis tool, but I think it is not accurate enough. is there a way in 3D analyst or Spatial analyst?
Thanks

In the case that roads overlap each other on different elevation, is there better way to still find the elevation for each street centerline segment (from node, and to node)? I have both contour line and DEM for the area. But the roads overlap a lot.


texts:
I am trying to create a Delaunay triangulation of a set of 3692 points with QGIS 1.8.0. The resulting set of polygons appears unfinished, leaving a small part of the original map uncovered by triangles and a small set of points not used in any triangle. Is there a limit of the number of input points or is there some other constraint?

Best regards from France


texts: QGIS 3 access "extract by location" processing tool with python

texts:
I am trying to use the "Extract by location" algorithm in a standalone script with QGIS 3.4.7 and Python 3.7.
The algorithm, however is not listed in available algorithms as listed with: processingRegistry().algorithms().
Since there have been a lot of changes in QGIS 3, and the docs are not complete, I can't figure out how to use it.
The algorithm exists in the app, since I can access it and use it from the QGIS app itself.
Of the algorithms that have anything with 'location' or 'select', I can only use:

qgis:joinbylocationsummary
qgis:joinattributesbylocation
qgis:randomselectionwithinsubsets
qgis:selectbyattribute
qgis:selectbyexpression
qgis:extractspecifiedvertices

So I guess the processing tool for "Extract by location" is not using a dedicated algorithm object anymore, and it is using some other tool. Can anyone suggest, how to check what method/class/object is used to run "Extract by location" in QGIS 3?

In QGIS 2.x the algorithm was used with
processing.runalg(...) but this is not used anymore, so all guides for running algorithms with python now do not apply.

My code starts with

qgs = QgsApplication([], False)
qgs.setPrefixPath(r"D:\QG\qgis-ltr", True)
qgs.initQgis()

sys.path.append(r'D:\QG\qgis-ltr\python\plugins')
from processing.core import Processing
Processing.Processing.initialize()

I use this code to list the algorithms:

qgs.processingRegistry().algorithms()

And this code to run one:

Processing.Processing.runAlgorithm('qgs:extractbylocation', <args>)

Which throws:

    raise QgsProcessingException(msg)
_core.QgsProcessingException: Error: Algorithm qgs:extractbylocation not found

as if the algorithm was not available.

The QGIS app allows to use it with no problem.
.


texts: Counting where multiple street names within polygon using MapInfo?

texts:
I have a dataset of address points with their associated street name, in addition to a tab file of polygons covering a catchment. I want to know if there is a way of selecting polygons which have address points with different street names within them.

I am guessing that there is a way of using the count function in MapInfo but unsure how this works with character fields.


texts: Publishing geoprocessing tools to ArcGIS for Server gives ERROR 001488?

texts:
I have been trying to publish a geoprocessing tool to our new ArcGIS for Server 10.1 with SP1. I configured the serves to not to copy data to server when publishing services.

The tool ran well and was shared as a service definition file successfully from ArcCatalog. It however always failed in publish server from ArcGIS Server Web Manager due to the following error:

  ERROR 001488: The service you are trying to publish includes data but
  the site has been configured such that data cannot be copied to the
  server. ERROR 001364: Validation of service definition failed. Failed
  to execute (Publish Service Definition). Failed.

Can anyone help with any hints?


texts: Computing total road length in fixed radius circle around a point

texts:
I have a dataset of a number of points which are plotted on QGIS using their latitude and longitude. Now, I have to find out the total length of roads surrounding each point within a circle of radius, say 100m.

This feature comprising of total road length surrounding each point in a buffer radius is further to be used for land-use regression.

Can someone suggest on how to obtain this feature?


texts:
I need to plot a map to show monthly SO2 emission point over a specific area. To do so, I have an Excel monthly dataset (derived from satellite data) contains a column for latitude, a column for longitude and a column for SO2 amount. However, when I add XY data on ArcGIS, there are lots of SO2 values for each lat/lon (because satellite scans the area every day, and I add data for 30 days at once).

The problem is that when I use an interpolation tool (kriging), I do not get the right result. How can I sort this problem out? The image I tagged is related to when I add XY data on the map.


texts: Projecting coordinates into different EPSGs gives inconsistencies

texts:
I am using 3 EPSGs: 4326 (WGS84), 3857 (Mercator) and 3035 (EU Copernicus dataset).

I am trying to easily convert coordinates from one system to any other. To do so, I wrote a Coordinate class with the following method, which serves to take the point's coordinates self.x and self.y given in self.epsg and project them into new_epsg:

   # Project point into given EPSG system
   def project(self, new_epsg):
      # Create transform from src epsg to dst epsg
      src = osr.SpatialReference()
      src.ImportFromEPSG(self.epsg)

      dst = osr.SpatialReference()
      dst.ImportFromEPSG(new_epsg)

      transform = osr.CoordinateTransformation(src, dst)

      # Create and project point
      wkt = 'POINT ({} {})'.format(self.x, self.y)
      point = ogr.CreateGeometryFromWkt(wkt)
      point.Transform(transform)

      self.epsg = new_epsg
      self.x = point.GetX()
      self.y = point.GetY()

      return self

However, when converting from an epsg to another the coordinates are always wrong.

Testing I've done:

I've tried a combination of two things:

flipping the coordinates in input i.e. wkt = 'POINT ({} {})'.format(self.y, self.x) instead of self.x, self.y
flipping the coordinates in output i.e. self.x = point.GetY() self.y = point.GetX()

To test my code, I've used a single point in the 3 ESPGs using https://epsg.io/. The coordinates are (x, y):

(1) EPSG 4326: (-31.22000, 39.443505) see here https://epsg.io/map#srs=4326&x=-31.22000&y=39.443505&z=12&layer=streets
(2) EPSG 3857: (-3475394.537268, 4785400.320095) see here https://epsg.io/map#srs=3857&x=-3475394.537268&y=4785400.320095&z=12&layer=streets
(3) EPSG 3035: (949571.84, 2769907.62) see here https://epsg.io/map#srs=3035&x=949571.84&y=2769907.62&z=12&layer=streets

Now here are the results:

not flipping anything: nothing works
flipping only input: only (1) to (2) and (3) to (2) work
flipping only output: only (2) to (1) and (2) to (3) work
flipping both: only (1) to (3) and (3) to (1) work

So overall by flipping coordinates around I can get all cases to work but this is very case specific and I would much rather have a general one-fits-all function.

More (runnable) code is available at https://pastebin.com/vzDkC3QM

I hope I didn't just miss something stupid.


texts:
I am able to download osm data with the following url request
"http://www.overpass-api.de/api/map?*bbox=7.505549,6.482084,7.511999,6.488254". It is generating a map.osm file, in that some of the roads are extending beyond the bbox limit and not exactly clipped to the bbox.

How do I sort this issue?


texts: Creating a Network Dataset in ModelBuilder?

texts:
I am working on modeling a project that requires me to create a network dataset in ArcCatalog. To do this, you simply right click on the roads shapefile and choose, "Create Network Dataset". You then just follow the wizard. I was wondering how that can be modeled in ModelBuilder.

There is a "Build Network Dataset" tool, but that doesn't seem to be what I'm looking for.


texts:
What is the best way to detect spatial outliers in R? With spatial outlier detection I mean:

  Discover the objects whose non-spatial attribute values are
  significantly different from the values of their spatial neighbors.

Maybe performing univariate Z algorithm, Moran scatterplot or Scatterplot approach like this paper describes: On Detecting Spatial Outliers?

Or performing Metric Learning like this paper describes: Contextual Spatial Outlier Detection with Metric Learning?

Other approaches are welcome!
I'm looking for a working implementation using R code.