texts: I like to know how can I color a xyz file based on the Z values(elevation) with QGIS, by python or graphicly


texts: Creating a time-series from an image-collection in GEE including null-values in case all pixels are masked

texts:
I want to create a time-series from an image-collection resulting in a list with the time-series values. It could be however, that some of the images in the collection are completely masked (i.e. have no valid pixels) in my Area-of-Interest. I create the time-series in GEE using this script (the first line is just to make a dummy image-collection, which has one fully masked image) (https://code.earthengine.google.com/9aa3e051390d72fac3fa783b009b08b4):

var img_col = ee.ImageCollection.fromImages([ee.Image(5),
                                             ee.Image(7).updateMask(ee.Image(0)),
                                             ee.Image(10)]);

var create_ts = function(img){

  var value = img.reduceRegion({reducer: ee.Reducer.mean(), geometry: aoi, scale: 10000});

  var ft = ee.Feature(null, value);

  return ft;
};

var ts = img_col.map(create_ts);

print(ts);
print(ts.aggregate_array('constant'));

The image-collection contains three images, the printed ts shows only two values however:

[5, 10]

but I would preferably like to have something like this:

[5, null, 10]

as output. If it's not possible to have "null" values in a list/array, then some kind of No-Data-Value would be acceptable as well, like for example:

[5, -9999, 10]


texts: How to deactivate zoom on the mouse wheel?

texts: I need to fixed scale in canvas in QGIS but I didn't find this function. If I fix scale it is impossible to make zoom in or zoom out into canvas. Only one move which layers loaded can be do are top/down and left and right.


texts: Exporting QGIS shp file to dxf-Which dxf format?

texts:
As I understand it, dxf comes in a range of format types. QGIS (v2.0 'Dufour') can export any layer to a dxf file. However I have a hydraulic modelling package which relies upon importing the dxf file as the R12 format.

This import does not seem to work. I was wondering which dxf format does QGIS export to? If it is not R12, this may be the reason why this is not working.


texts: Converting from British National Grid (27700) to WGS84 (4326) to Web Mercator (102100)

texts:
I'm developing a web application that uses ArcGIS to display a map using JavaScript.

We are holding all our data in MS SQL Server 2008 as Easting/Northing against the British National Grid. Unfortunately ESRI needs the data plotted to Web Mercator.

So I've ended up using Proj.Net for the British National Grid (27700) to WGS84 (4326) conversion

Then for the WGS84 (4326) to Web Mercator (102100) conversion I used the code in this article.

This gives me usable data within ESRI, however plotted points aren't quite in the right place. It seem that I'm suffering from the shift because the datum hasn't been changed. This post suggests that Proj.Net can correct the datum shift by specifying the TOWGS84 parameter, but I don't understand what this means.

I'd prefer a C# library or code that got from British National Grid to Web Mercator in one quick, simple conversion. However if I need to jump through hoops, I'll do that, but could do with some extra guidance


texts:
I am working with a wildlife-related feature class where I have different fields for common name (i.e. Northern Goshawk), scientific name (Accipiter gentilis), and species code (NOGO). I have domains for all of these which limits acceptable answers to only the species that occur in the area studied. However, as I'm essentially adding the same information 3 times per data point, any species code I enter will only match up with one common name and one scientific name. Is there a way to link these fields so when I enter one field information (say, NOGO) the other two auto-fill?

I'm quite new to geodatabases.


texts: I have more than 100 QLR files, QGIS does not allow me to select multiple and, adding them one by one takes a long time. How can I add them all at once?


texts: Tool to set the maximum shape area of drawn polygons [ArcGIS 10]

texts:
Is there any tool or method to define beforehand the maximum area of each drawn polygon class ? Polygons are going to be distributed according to my idea so I need to have control over the process. I can imagine it the way that while drawing polygons I am not able to exceed the previously defined threshold value.

I started doing it manually (had to create patches of tree species of exact total size in ha) and was checking the total size in a table of attributes but I believe there is a way to do it in more comfortable and precise way.


texts: How can I get latitude and longitude from FIPS codes?

texts: I have full FIPS codes down to the block, looking for a way to convert the FIPS codes into centroid latitude/longitude for each block. Basically this but backwards (and I'd like blocks, not just tracks): API to reverse geocode latitude, longitude to census tract?


texts: GDAL: *.nc rasters, opened as "netCDF/Network Common Data Format" under Linux and as " HDF5/Hierarchical Data Format Release 5" under Windows

texts:
I have a set of rasters that are in the "*.nc" format and I noticed that gdalinfo (and subsequently other GIS packages based on GDAL like QGIS) is opening them with two different drivers.

Linux:

gdalinfo A20133052013312.L3m_8D_CHL_chl_ocx_9km.nc
Driver: netCDF/Network Common Data Format

Windows:

gdalinfo.exe A20133052013312.L3m_8D_CHL_chl_ocx_9km.nc
Driver: HDF5/Hierarchical Data Format Release 5

The problem is that under Linux I get the proper/expected range of values, while under Windows the range not only is not the same, but is completely off.


texts:
I am trying to use Timemanager with my csv files which i loaded on my map as shp layer. those contain time staps which i tried in hh:mm:ss.fff hh:mm:ss:fff hh:mm:ss,fff and also each with 1 f - in the exvel csv and in timemanager and layer properties.
Time is running but layers arent displayed.
I then loaded the data without milliseconds hh:mm:ss and it was working after i formated the excel table to that time stamp.

However, i still cant get milliseconds to be run with timemanager.
Can anybody help?


texts: Writing automated tests for QGIS plugins

texts:
I'm looking for advice on writing automatic tests for QGIS plugins written in Python.
I have written tests for Python scripts in the past using PyUnit (the unittest module), but have never done so for an application with a GUI. I've found a page describing how to use PyQt4.QTest to do unit tests on Qt widgets (http://www.voom.net/pyqt-qtest-example), but I'm struggling to see how I can use this with a widget that has been designed to run from within QGIS.
The section on "Testing" in the PyQGIS documentation is notably absent.
What I have so far is:

Keep the actual processing of data in isolated modules or functions, and write unit tests for those;
Perform basic testing of the UI using QTest;
Pray that it all holds together when using the plugin from within QGIS.

Is there a better way?


texts: How long should it take to export data to your drive?

texts:
I'm new to GEE and I'm using their test script to send an image to my drive:

// Load a landsat image and select three bands.
var landsat = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_123032_20140515')
    .select(['B4', 'B3', 'B2']);

// Create a geometry representing an export region.
var geometry = ee.Geometry.Rectangle([116.2621, 39.8412, 116.4849, 40.01236]);

// Export the image, specifying scale and region.
Export.image.toDrive({
    image: landsat,
    description: 'imageToDriveExample',
    scale: 30,
    region: geometry
});

My first attempt ran for 20 minutes with never fully downloading.  The second attempt crashed after about 20 seconds and kicked out a earth engine exceeded memory capacity error.  I'm currently running it again, but so far it's almost 20 minutes again.

Can some one tell me how long on average it should take, and if I'm doing something wrong?


texts: Relationship Class Table displays empty table

texts:
I'm trying to link 2 tables together using relationship classes (pre-defined within this datamodel). I've uploaded a picture which might be clearer, but here is the process I'm trying:

Table named Site, with a field EventID.

Table named Anode, with a field SiteEventID.

EventID and SiteEventID are the same (same values).

1-M relationship class (simple) / Origin name = Site, Primary key = EventID, Foreign = SiteEventID / Destination name = Anode / Labels forward = Anode, Backward = Site. Relationship name = SiteAnode

When I open either Site or Anode using SiteAnode, the table that opens is completely empty.

The table remains empty until I remove it from ArcMap, and reload it from the gdb.

Has anyone experienced this, or am I missing something?
Something else to mention - if I create a feature from Anode, and perform a Join with the same key it works fine.


texts: Leafletjs : combine location data in geoJSON/KML/GPX with other information in CSV through key field

texts:
I've been able to load CSV having lat-long co-ordinates into leaflet using this useful solution. Looking for ways to cross-link co-ordinates from a different source and data from another. Describing the use case:
Polygons representing wards in a city are loaded through geoJSON or KML or GPX file, each object having a unique name like "Ward 01". This file has the polygon co-ordinates of each unit, and its name, but not much else. Now, there is another file, a CSV, like this:

name,detail1,detail2,detail3,detail4,detail5
Ward 01,23,45,10,54,30
Ward 02,78,98,14,20,53

..and THIS has the details that we want to see in a popup when the user clicks on a polygon, or mouses over and we display on an info panel, etc etc.
How does one combine the two? So far I've seen solutions where the location data and the details are all in the same file. Real-world cases aren't that simple ;)

The "name" field is the unique non-repeating connector here between the geoJSON/KML/GPX and the CSV.

An alternative is to merge this data sitting in the CSV into the file having the spatial information. Is there an automated way to do it? One hack is to convert the mapping data to a tabular row, merge the data in and convert it all back to geoJSON (or do the same with KML/GPX), but I'm a little vary of taking that complicated route, and it will hamper the occasional updates to data (by non-coders) that this use case requires. It would be a lot more efficient to just pull the details from the CSV for each object in the HTML itself.

Addendum : Yes, using MySQL DB instead of CSV is also acceptable though not preferrable (serverless, portable operation desirable). But please share how to go about it, because I still don't get how to pull it in.

Edit: One major hurdle I'm finding here is the asynchronous loading of the data. Most solutions I see have the entire code loaded into a file-loading script, to be called after the data has finished loading. But here, there's a main program with multiple overlay layers, and a CSV each to load for each layer.

Edit: Here's a link to the project, see source code:
http://nikhilsheth.techydudes.net/files/pune-pb.html


texts: Find line segment that is not between other lines with QGIS

texts:
As shown in below image, I have a project with rivers(blue) which are bordered by fences(green). I have to identify the places that are not guarded by fences. The fences can be further away from the water, anything over 20 meters away does not count as guarding the water anymore.

Im using QGIS, also open for Python solutions (Fiona/Shapely/Geopandas/..)

My best try:

Create 4 buffers on the fences: North, East, South, West, using MMQGIS
Find overlap of North & South, West & East.
Find difference between rivers and fences

However this fails because MMQGIS can appearently no longer create buffers without errors.


texts: Weights assignment based on intersection

texts:
I'm trying to assign weights for the grids in postcode based on intersection with urban layer for each postcode (sum of weights must be 1). To get the count of intersection between grid and urban layer I have created a temporary table (vrg.urban_nonurban_count_new ) which gives me postcodeid,urbancells count, non urban cells count.

For example in the screenshot below:
 There are 5 grids in this postcode(postcode containing centroids). When I execute the above query, I get the below result

Below is my query:

SELECT
   ST_Union(x.geom), x.vrgid, x.weights
FROM
     (SELECT
           postcodes.id as postcodeid,
           fishnet.gid as vrgid,
           fishnet.geom as geom,
           CASE WHEN st_intersects(centroids.geom, urban.geom) and
                 counts.nonurbancells!=0
            THEN 0.95::numeric/counts.urbancells
            WHEN st_intersects(centroids.geom, urban.geom) and
                 counts.nonurbancells = 0
            THEN 1.00::numeric/counts.urbancells
            WHEN st_disjoint(centroids.geom, urban.geom) and
                 counts.nonurbancells!=0
            THEN 0.05::numeric/counts.nonurbancells
            END  AS Weights
FROM vrg.urban_nonurban_count_new as counts
LEFT OUTER JOIN vrg.gfk_2016_id_5_digit_pcd_areas2013_projected as postcodes
            ON postcodes.id=counts.postid
LEFT OUTER JOIN vrg.rdsid_86_quadgrid_centroids as centroids
            ON  st_contains (postcodes.geom,centroids.geom)
LEFT OUTER JOIN vrg.rdsid86_quadgrid AS fishnet
            ON  fishnet.gid=centroids.gid
LEFT OUTER JOIN vrg.rdsid86_katrisk_poly_projected as urban
            ON st_intersects(urban.geom,fishnet.geom)
WHERE postcodes.id='42395') as x
GROUP BY x.vrgid, x.weights

Here I'm supposed to get 0.0125 for 4 cells and 0.95 for 1 cell. But, 2 weights have been assigned to single cell i.e 0.95 and 0.0125 for vrgid 7192.
I'm not able to figure out the problem.


texts: Writing hidden information to PDF file from ArcMap?

texts:
I'm usually creating maps for (different aspects of) reports of large infrastructure projects. Layers and attributes tend to change back and forth as time goes by, and I'm forced to create new versions of these maps regularly.

The problem is that in each project I have some 50-60 mxd files, some with very small changes. What makes it even more complicated, all of the produced maps should have the same date (and an overall very strict layout) according to that specific project's end date.  Even with a very thought-through naming convention, it can sometimes be a real pain to find out which mxd belongs to a certain map that I am tasked to update.

Is there a way to automatically write the path of the mxd (and preferably the date of export) to a pdf and/or a png file somewhere in the metadata when exporting from ArcMap?

The information should not be visible when printed, nor could the file name be used for this (since it will change multiple times), but should also not be too complicated to find when looking for it. The process of exporting a pdf/png in this way should be reasonably fast, since I'll be doing it several times a day.

Does such metadata exist within pdf files and/or png files, and how can I access it?


texts: Is there any difference in measurements when using the draw toolbar vs the editor toolbar?

texts:
Say I have a cornfield that I wish to measure acreage on.

If I utilize the draw toolbar, create a polygon around the field, export to shapefile, add area field in attribute table, and calculate geometry for area in acres, will this be accurate?

I have been using the editor toolbar to edit and cut polygons to shapes that fit my needs and then calculating geometry for acreage on those, but I am looking for an easier way and wondering if this process of using the draw toolbar and creating polygons was synonymous with using the edit toolbar.


texts: Find nodes inside area (like forest) with OverpassAPI

texts:
I'm still in the processing of learning how to use the OverpassAPI QL syntax in regards to areas.

With the following query I get the forrest areas:

[out:json]
[bbox:{{bbox}}];

(
  way["landuse"="forest"];
  relation["landuse"="forest"];
  way["natural"="wood"];
  relation["natural"="wood"];
);

out geom;

With this query I get the nodes I'm looking for inside the area queried before:

[out:json]
[bbox:{{bbox}}];

node[natural=spring];
out;

When I combine those to I don't get any output (in Overpass-Turbo):

[out:json]
[bbox:{{bbox}}];

(
  way["landuse"="forest"];
  relation["landuse"="forest"];
  way["natural"="wood"];
  relation["natural"="wood"];
);
map_to_area;
node(area)[natural=spring];

out;


texts: Reading through contents of layer file using ArcPy?

texts:
I have a collection of rasters (all in separate geodatabases) that I am trying to put into a MOSAIC using a python script. I have a bunch of these to make so I was thinking that the easiest way to do it would be to create a .lyr file containing all the rasters I want to add and then have a script read through that layer file and add the rasters using AddRastersToMosaicDataset_management. The problem is that I haven't been able to read through the contents of a layer file. This line seems to break my code: layerfile=arcpy.mapping.MapDocument(r"S:\path_to_layer.lyr"). I try to catch the error message with except: print arcpy.GetMessages() but nothing prints.

I thought I'd try the other approach and use an mxd instead of a layer file. I created an mxd with just the rasters I wanted to add, but I couldn't read the path of the raster (shapefiles did work though). Here is my code:

import arcpy
mxd=arcpy.mapping.MapDocument(r"S:\path_to_mxd.mxd")
dataframes=arcpy.mapping.ListDataFrames(mxd)
for dataframe in dataframes:
    for lyr in arcpy.mapping.ListLayers(mxd, "", dataframe):
        desc = arcpy.Describe(lyr)
        print desc.catalogPath
        #print(arcpy.Describe(lyr).catalogPath)

This will print a path to a shapefile, but gives IOError: "Raster_Name" does not exist when the layer is a raster.

At the moment my only working solution is to manually generate a text file with the path to each of the rasters and read through that line by line. I don't want to have to generate a text file for every MOSAIC I need to create. I feel like I've got 4 possible paths I could go down and I'm getting stuck on each one. Can anyone help me either;

read through a .lyr file or an mxd and extract the paths of the raster datasets; or
convert a .lyr or an mxd to a text file listing the paths of the datasets contained?


texts: Adding layer names to combobox in QGIS plugin - layer names appear more than once

texts:
I am creating an OGIS plugin, but I have a problem. I want the layers in the problem layer panel to automatically appear in the ComboBox. I wrote the necessary code for this, but when I open the plugin every new layer, each layer in the layer panel appears twice, how can I fix this?
Code written:

Layers = self.iface.legendInterface().Layers()
Layer_list = []
For layer in layers:
      Layer_list.append (layer.name())
      self.dlg.comboBox.addItems(layer_list)


texts: Best way to manage large number of GeoTIFF rasters?

texts:
I have around 2-3000 1-bit GeoTIFFs of size 11000x9400. Color 0 should be transparent, and QGIS gets it wrong in the first place, so I used a .vrt file to set all the parameters in one single place.

Problem is, even if I set the layer to be visible only at a small scale, the project takes a lot of time to load in QGIS.

Which is the best way of dealing with all these files? I tried using grass, but I don't see if it can really be faster than QGIS.


texts: How to repair broken Mosaic Dataset data sources?

texts:
I have a question. My computer has been replaced with a new Xeon Workstation. Before that, I had a temporary workstation and had ArcGIS 10 on it. While I was waiting for the new workstation to be replaced, I was working on mosaic dataset and worked on several clipped rasters to mosaic them. When I got my new workstation back, I had to backup all of my datasets to external hard drive and then load them to the new computer.
So when I bring the mosaic dataset up to my ArcGIS 10 on my new computer, all I see is black! Did I do something wrong when I made the backup or should I run some fixes?

Originally, what I did was that I had two NAIPs loaded from the ArcGIS server added them to my ArcMap and had a layer called mask to clipped the raster within the mask. I use them to mosaic them to  created  one mosaic file saved in Mosaic Dataset.

HELP !! I hope I do not have to re-do all the work again. Thanks for your feedback.

Thanks !


texts: I want to create a tool that will search a single field, in a single layer and zoom to that location?

texts:
I have a feature class of valves and another of fire hydrants. Each has a unique identifier field. I want to create a tool that has one dialogue that allows input of a identifier then zooms to that location. I have a solid understanding of coding and GIS. But don't even know where to start with this.

Any clues about where to start with this?


texts:
I'm using a new install of QGIS 2.18.2 on an iMac running OSX 10.10.5.

I create a new project.  I use the OpenLayers plugin to show an OSM background of a part of Mayfair in London.  I now use the Georeferencer to add an 9Mb .png (copy of a web map) and import three points. All the SRS/CRS are set to EPSG:3857. Transformation settings are linear.

This Dropbox folder contains a copy of

the 9Mb map I'm trying to transform
the points I have previously saved
a screen copy of the transformation settings
a screen copy of the QGIS window showing the GCP table and the two maps

At this point I run the transform and I get the message:
Transform failed: Failed to compute GCP transform: Transform is not solvable.

I have no idea of what is going wrong here.  I have managed to use exactly the same process (with a different but very similar .png) and it has worked, but an attempt to add a second file usually seems to fail.   Trying to add the map that did work to this project also failed in the same way.  Many other images/points have failed as well.   I've also reduced the PNG to a 1.6Mb jpg and this also fails in the same way.


texts: Standard deviation and average over 1000 images with random null values

texts:
I have a large collection of raster files that I would like to calculate the standard deviation over. The first issue is that ArcGIS cannot process over 1000 rasters via cell statistics, so I decided to go another route by summing the squared values of all the rasters, and then getting the variance and standard deviation. The problem is that each raster may have different areas with values, and some areas with null.

If I were to try to script an iterative function that would sum the rasters, the output would be completely null since with a large quantities of these maps would eventually have null values everywhere.

I then tried using another program, Spirits (http://spirits.jrc.ec.europa.eu/overview/about/). While this software was able to provide an image like this:

The problem with this output is that the null values were converted to 0.

I am looking for a method to process a large time series of these images that ignore the null values for their calculations, and can process over 1000 images at once. Do you guys have any recommendations?


texts:
I am using a custom map with the TMS in OpenLayers.  I know my high/low latitude and longitude points.  In other words I know the exact bounds of my map image in Lat/Long.

Is there a way to use lat/long in the bounds?  Or do I have to convert to Mercator?  My bounds are pretty small - it only covers maybe 10 square miles and then I want to be able to draw circles on the map at specific lat/long locations - like draw a circle to highlight a gas station.

I am just a bit confused about the Bounds portion.  Here is my code this far - use they Mercator bounds just to get my image to display:

function loadmap() {
        var get_my_url = function (bounds) {
        var res = this.map.getResolution();
        var x = Math.round((bounds.left - this.maxExtent.left) / (res * this.tileSize.w));
        var y = Math.round((this.maxExtent.top - bounds.top) / (res * this.tileSize.h));
        var z = this.map.getZoom();

        var path = z + "_" + x + "_" + y + "." + this.type;
        var url = this.url;
        if (url instanceof Array) {
            url = this.selectUrl(path, url);
        }
        return url + path;
    }

    var options = {
        maxExtent: new OpenLayers.Bounds(-20037508.3427892, -20037508.3427892, 20037508.3427892, 20037508.3427892),
        numZoomLevels: 5,
        units: 'm',
        projection: "EPSG:900913",
        displayProjection: new OpenLayers.Projection("EPSG:4326")
    };

    var map = new OpenLayers.Map(
    'wissmap',
    options
);

    var layer = new OpenLayers.Layer.TMS(
    'Aerial',
    '/custom/maps/',
    {
        type: 'jpg',
        getURL: get_my_url
    }
);

    map.addLayer(layer);
    map.setCenter([2048, 2048], 2);
}


texts: ERROR: parse error - invalid geometry (Postgis)

texts:
Initially I was storing my column as path and then following this question I changed it to text.

I'm trying to convert my column from text to geometry type by running this command:

ALTER TABLE area_grid_01 ALTER COLUMN coordinates TYPE Geometry(point, 4326) USING ST_SetSRID(coordinates::Geometry, 4326);

And facing this error:

ERROR:  parse error - invalid geometry
HINT:  "[(" <-- parse error at position 2 within geometry
SQL state: XX000

Following this question, I changed my data formatting and column looks like this:

The reason I am converting it into geometry in the first place is to be able to perform postgis functions on data.

I have explored many other ways trying to figure out a way around it. Seemingly, its not working yet!


texts: How to disable schema locking on all services on ArcGIS Server?

texts:
I have a situation in which all services published to an ArcGIS Server (10.6.1) were published with schema locking enabled. I would like a way to disable schema locking on all services without having to go into each one and set the schemalockingEnabled parameter to false.

The source datasets are mostly from an enterprise geodatabase (SQL Server 2017).

I'm looking for a python solution, possibly using the ArcGIS API for Python. Unless there's a non-programmatic way to accomplish this.

Is this possible?


texts: Implementing Equal Earth Projection for ArcGIS Desktop?

texts: I am trying to find a means to implement the new Equal Earth projection (http://shadedrelief.com/ee_proj/) into ArcGIS Desktop (for my contractors to use in our products) - has anyone made progress on this yet and where could I access this?


texts:
Over the past several days I've read as much as possible regarding the OpenLayers API (within a JavaScript/jQuery environment), yet I'm still having trouble in getting my head around what I believe to be basic features.

I can successfully place maps on my page, set bounds, etc, etc, but I'm having an incredibly hard time in the following:

My requriements are extremely simple, I want a 'file' (KML, JS, I don't know) of each country and it's borders so that I can simply (as an example) select a country from a dropdown list and OpenLayers highlight that country!

I did think it would be a simple task, have gone through all 220 examples of the OpenLayers website and can't see anything remotely close.

I've been here: http://thematicmapping.org/downloads/world_borders.php

Downloaded the WorldData which I thought would have the vector co-ordinates for the worlds countries borders, it looks to be a database file but haven't a clue what to do with this!

If someone can break down the basics, I'd be incredibly grateful and better still show a basic example of what I'm after, I should be able to learn from it.

Many thanks,


texts:
My application is throwing the following error:

DatabaseError: GetProj4StringSPI: Cannot find SRID (3857) in spatial_ref_sys

I'm guessing (as stated above) that I need to add the SRID 3857 to the spatial_ref_sys.
Looking at my notes from a previous project I've accomplished this in the past by doing the following:

Adding an entry for that SRID to the epsg file, then

Running the following from python shell

from django.contrib.gis.utils import add_postgis_srs

add_postgis_srs(3857)

My problem then becomes, I don't know what the epsg entry should be for 3857.

Am I approaching this problem correctly?

If so, does anyone know what the epsg entry for 3857 should be?

Is it just the same as 900913 (also spherical mercator)?
In other words, could I just replace 900913 with 3857 in the following?

<900913> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs


texts: Can't access GeoServer after CORS filter

texts:
To enable CORS filter I added the following in web.xml file in webapps/Geoserver/WEB-INF

<filter>
 <filter-name>cross-origin</filter-name>
 <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>

...

<filter-mapping>
    <filter-name>cross-origin</filter-name>
    <url-pattern>/*</url-pattern>
</filter-mapping>

After adding above, I restarted GeoServer. Now, when I access GeoServer, I get following error:

HTTP ERROR: 503

Problem accessing /geoserver/. Reason:

Service Unavailable


texts:
Using FME Desktop 2017.1.0.0 (20170731 - Build 17539 - WIN64) reading a Excel xlsx file on Windows fails and returns the following error:

  ...|ERROR |Excel Reader: Could not open
  file '...\input.xlsx', file is too large. Try switching to 64-bit FME,
  or exporting the Excel file to CSV. If the file is exported to CSV,
  formulas and formatting will be lost

The file is 56KB, has 2 sheets with a couple of rows and columns. The data is imported with a Query from 2 csv files, what might cause the problem.

Converting to CSV is not an option as this is part of an repetitive automatized process.

Copy pasting the data to a new excel document solved my problem but I would like to keep the source data original. Does anyone have the same experience and found a solution?


texts: adding helicopter video files to Google Earth or QGIS

texts:
Is it possible to add a helicopter's video files to Google Earth or QGIS ?
We have a flight path points derived from the xml file associated with the .mpg file which I envision being used to pilot the video.
The .mpg files range in size 100mb to 2gb.


texts: need to show multiple symbols for one point

texts: I am working on a state map of parks with campsites. I want to show symbols for handicap access, tent camping, backpack camping, winter camping, rv camping. I have used a symbol for handicap but it shows up multiple times at some parks. Also if there is more than one type of camping it is over shadowed by handicap. any ideas? thanks


texts: How to do separability analysis in eCognition?

texts: Can I use eCognition to check which bands are more helpful to separate individual classes and to get output result? If so, how can I do this?


texts:
I am new to CartoDB and NOT a programmer. I am trying to make a simple torque map that changes by year and I do not want the slider text to include or display the month, just the year.

Is this possible to do with little to no programming?.


texts: Tool to colorize Pointcloud from 360° Laserscann

texts:
I have a Pointcloud with panoramic images that I want to colorize. What tools could I use for this. I looked into PDAL but from what I read, PDAL needs orthographic pictures.
Is there another tool worth looking into?


texts: I am currently using Geocortex Essentials REST elements and I am in the process of adding a new map service. I've entered in my url to the map service's in ArcGIS Server, but can't see my newly created map service in the dropdown list. Any ideas or hints as to why my map service isn't appearing?


texts: Can I export a polygon shapefile to a Garmin device?

texts:
I want to export the boundaries of a parcel to my GarminMap 60CSx. I downloaded a parcel map from the county, and I want to use it to walk the perimeter.

Does the .gpx format support polygons? Or is there some other way?

I'm on Mac, using QGIS 1.7.4.


texts:
I'm using the qgis2web plugin to display a Landsat image and some agricultural fields. Both raster and vector data are in EPSG:4326.
In the QGIS canvas, the vector and raster line up (left side of the attached image) but in the qgis2wweb dialog there is a shift in the raster (indicated by red line). I know it's the raster that's shifted because the vector lines up with OSM background. Has anyone seen this?


texts: Why are poles missing from GIS topographic data?

texts:
Why are poles data missing from most topographic resources ? Data from 60⁰ to 90⁰ lat frequently missing from topographic GIS data. Example with SRTM4.1:

Common sense or deduction lead me to think that either:

Orbits: the orbit does not go to poles. Which is absurd for a topographic mission.
Data : the data on poles is messed up. But come on, in 2010's, space agencies computer scientists are not able to clean this up ?

I asked myself this for years. Yet, I have no real, solid answer up to now.

Just for fun, if the there are relevant GIS resources for these 2 areas, please share.

Edit :
- University of Minesota: REMA (Reference Elevation Model of Antarctica). Specifics: topography including snow layer(s), 1px= 8m, 60-88° South (Antartica),. Files: 2x2m;8x8m;...,. —♣ Comment(s): Full data is 43TB large!


texts: Algorithm or library for buffer zone for lat-lon (wgs) coordinates in global scale

texts:
Can anyone point to generic algorithm or library for buffer zones if the source data is in WGS84 degrees? The buffer width should be in meters, as buffer in degrees it generally does not make much sense.

I'm aware of popular geometry libraries like JTS, PostGIS, ESRI etc, but none of these seem to work for lat/long case. I can of course convert latlon to some cartesian projection and then use these libraries, but this is usable only in small scale, limited specific area and/or with sacrificing precision. Also for some cases you may "convert" meter distance to degree distance, but this is ugly and heavily dependent on latitude. Some libs have first algorithms for proper lat-long, e.g. PostGIS can do Haversine distance, but no Buffers.

So basically I would need a Buffer polygon which would be created using Haversine distance from given line.

Now if you really do not know such algorithm/tool, which would work without limitations (globally, with any coordinates and buffer size, giving still high mathematical precision), then you can suggest "second class" alternatives with limited input area and buffer sizes, keeping precision to <1m scale and assuming that the area is somewhere in inhabited regions (i.e. not in Antarctica or ocean), and total size is below say 500km. For most cases it would do the job also.

p.s. tools and code samples in Java would be preferred. Specific target would be a mobile app (Android) which is used for technical underground infrastructure checks (digging safety) in cities anywhere in the world. It requires high precision regardless of latitude, but in sane range (<80 deg).


texts: Error when register database with ArcGIS server

texts:
I am using ArcGIS 10.1 (x32) and I have successfully connected to PostgreSQL(x64) as geodatabase.

When I try to register database with ArcGIS Server 10.1 (x64), I followed the instruction on the Esri' help. That is:

Download the x64 library for ArcGIS Server and copy to bin directory in install file.
Restart the ArcGIS Service
GRANT the database privilege to ArcGIS Server account

As I choose the connected sde database as register database, the error occurred:

The connection property set was missing a required property or the property value was unrecognized. Cannot connect to database because the database client software failed to load. Be sure the database client software is installed and configured correctly.

I am wondering anything else I should do to make this work?


texts: QGIS howto create raster from map canvas easily

texts:
i have question ho to create raster from QGIS map canvas (what i seen in my qgis) easily. Now i have used print composer and after some luck i get whole dataset into map item, after that i throw lucky guess for dpi and i get something.

What i would like to do is : make map look pretty in qgis, select cell resolution and get what i see in georeferenced raster file (geotiff, png etc..) from whole dataset. is this possible or do i just use print composer wrong ?


texts: Shapefile layers loading as the wrong geometry type

texts: I have a QGIS project shared with a few users and one user, in particular, is having an issue that when she opens the project the layers (shapefiles) are loaded incorrectly - they've lost their styles and the geometry is incorrect - lines appear as points.  Other users have no problems when opening the project; layers and styles are displayed properly.  We can sort it out by using the "change data source" and re-map the layers, but obviously, we shouldn't need to do this.  Are there any known issues that would cause this?


texts: R and sf - How to know this exists, data(us_states), without first using, (requireNamespace("sf", quietly = TRUE))

texts:
I am not new to using R, but I am newer to using the sf package, or any other similar geospatial packages.
I can easily explore and subset the world dataset and then plot the map of the US -  which includes Alaska and Hawaii. This data is easily explored from the top down using the pipe tool and subsetting.  But I actually only want to plot the map of the continental US.  Using the tools I normally use, I am unable to find an easy way of creating this plot.
A web search returned this website and showed the following code used to create the plot with the wanted feature:
if (requireNamespace("sf", quietly = TRUE)) {
  library(sf)
  data(us_states)

  plot(us_states["REGION"])
}

states49 <- (us_states["NAME"])
states49

What is this?  I've been using R for sometime and have never come across any code similar to this before.  I've never seen or used requireNamespace before. I never even knew this function existed.   As I mentioned earlier, it's easy to subset into the us_states using those tools I commonly use.
Using my regular tools, I never knew this data existed.  I would never have found this data without having used the code sample found from this website. Even the R official documentation says, "It's rare to load a package explicitly, but it can be done with:
requireNamespace()
loadNamespace() "

Link to quoted R documentation: http://r-pkgs.had.co.nz/namespace.html
So, to the point, how is someone supposed to know this data exists, when it might exist, when it doesn't exist, etc..  When would I suspect data like this might exist?  And when to find such data when such data might or might not exist?


texts: Clip bathymetry raster file from vector mask layer its not working in QGIS 3.4.3

texts:
I'm trying to extract bathymetry of Mediterranean Sea raster GeoTIFF file (downloaded from https://www.bodc.ac.uk/data/hosted_data_systems/gebco_gridded_bathymetry_data/gebco_30_second_grid/) from bathymetry (vector mask layer) downloaded from Natural Earth data.

After I selected bathymetry depth 0 at query builder of the vector layer, used Raster > Extraction > Clip raster by mask layer, but the area of elevation its grayed and not deleted.

I tried save as both layers raster and vector with same CRS WGS 84 and still the same result.
After save as the raster GeoTIFF (GEBCO2014_-12.3),(Gebco 4326 GeoTIFF) in CRS WGS 84 used again Raster > Extraction > Clip raster by mask layer.

How can i resolve this ?


texts: Difference between gdal versions in reading raster on two machines?

texts:
I am using raster data from Afripop.org and using the R package raster (raster_2.0-41) to do all my computation.  Everything has been working out just fine until I tried sending some jobs to a remote linux machine.  When I read in the raster on one machine the values are correct but on the other there are no NAs and the population (i.e. sum of all the values) WAY off.

On the linux machine I get:

> civpop <- raster("data/apciv10v3.flt")

> civpop[1:10]
[1]   0  60  28 198   0  60  28 198   0  60

> civpop
class       : RasterLayer
dimensions  : 7656, 7327, 56095512  (nrow, ncol, ncell)
resolution  : 0.0008333, 0.0008333  (x, y)
extent      : -8.599802, -2.494213, 4.357029, 10.73677  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : /home/andrewa/D4D/trunk/data/apciv10v3.flt
names       : apciv10v3
values      : 0, 255  (min, max)

One the other machine I (correctly) get:

R> civpop <- raster("data/apciv10v3.flt")
R> civpop[1:10]
 [1] NA NA NA NA NA NA NA NA NA NA
R> civpop
class       : RasterLayer
dimensions  : 7656, 7327, 56095512  (nrow, ncol, ncell)
resolution  : 0.0008333, 0.0008333  (x, y)
extent      : -8.6, -2.494, 4.357, 10.74  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : /Users/andrewa/Documents/D4D/trunk/data/apciv10v3.flt
names       : apciv10v3

If I use just rgdal I have the same issues so I guess it isn't an issue with raster. Any ideas?

Incorrect Machine (using GDAL 1.8.0, released 2011/01/12):

> civpop.gdal = readGDAL("data/apciv10v3.flt")
data/apciv10v3.flt has GDAL driver EHdr
and has 7656 rows and 7327 columns
> civpop.gdal@data[1:10,1]
[1]   0  60  28 198   0  60  28 198   0  60

Correct machine (using GDAL 1.8.1, released 2011/07/09):

R> civpop.gdal = readGDAL("data/apciv10v3.flt")
data/apciv10v3.flt has GDAL driver EHdr
and has 7656 rows and 7327 columns
R> civpop.gdal@data[1:10,1]
[1] NA NA NA NA NA NA NA NA NA NA`


texts: MapInfo tab files with LargeInt can't be converted to QGIS or any other format

texts:
I am trying to convert a MapInfo tab file with ogr2ogr
I did this:
ogr2ogr -f "ESRI Shapefile" testing.shp TEST.tab

And get this error;
ERROR 3: Failed to parse field definition at line 8 in file TEST.tab

    FAILURE:

    Unable to open datasource `TEST.tab' with the following drivers
    oid LargeInt ;

QGIS gives me this error when adding a vector layer:
Invalid Data Source: C:\Temp\Struer\TEST.TAB is not a valid or recognized data source.

I tried this:
Change the LargeInt to Integer in MapInfo.

Then it works fine.
So it seems like there is a problem with the LargeInt type.
Is it GDAL that can't handle this or does anyone knows a workaround?


texts: PyQGIS/QGIS 2.18: Filter layers by name in QgsMapLayerComboBox

texts:
I need a QgsMapLayerComboBox to show only layers whose name contains a specific word.
I have no way to know the names or geometry of the layers beforehand - I can only filter them based on whether or not they contain this word.

I thought about using setFilters(), but it looks like it only allows to filter by geometry or layer type, not by layer name. I also read about setExceptedLayerList() - but as I mentionned earlier, I don't know the exact name of the layers I don't want to show in the QgsMapLayerComboBox.

So, is there a way for me to use setExceptedLayerList() with an expression referring to layers whose name does not contain the word I need? Or another method that would allow me to filter the layers in my QgsMapLayerComboBox based on a word in their names?


texts: Zooming to point boundary using OpenLayers?

texts:
I'm wondering what the best strategy to position a map that has two points, using OpenLayers is. At the moment, my approach is to calculate the mid-point of the two points, and to zoom the map to that mid-point using the following approach:

var centerPoint = new OpenLayers.LonLat(lon, lat);
map.setCenter(centerPoint, zoomLevel);

However, this is not ideal as I don't know what the zoomLevel should be, so sometimes although the map is centered, the points are not visible.

Is there a better approach to dynamically position the map?


texts: Authorize app but also sign in user into ArcGIS online

texts:
I request authorization for my app through OAuth as described here:
https://developers.arcgis.com/documentation/core-concepts/security-and-authentication/browser-based-user-logins/

The sign in dialog appears with something like:

"App Name wants to access your ArcGIS Online account information".

But the sign in just gives access to the app. Confusingly it does not sign in the user which you can check by going in the same window to www.arcgis.com where you get sign in dialogue again.

Is there a way to actually sign in user into ArcGIS online? If I compare it with giving access by Google/Facebook etc, they give not only permission to app but also realy sign in the user.


texts: Failing to execute Quick Import of GML files to feature class database using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.6.1.

I am attempting to import GML files from the TOP10NL datasets (the historical data is only available as a package of GML files). In ArcCatalog, I selected "Dutch TOP10 GML" as the format and left the coordinate system as unknown. I selected all 14 GML files, and received the following message at the end of running the tool. The error, immediately following the message "translation was successful" and "Done import" was the message " Failed to execute (QuickImport).

I am using the Data Interoperability extension (authorized extension), and I have turned off background processing (as it is buggy for me).

[

So, the tool ran in entirety, but then failed. I'm not sure how to best address this, and I am really in need of this dataset.I don't know if it has anything to do with the size of the files or the number of features, I suspect this may be the case but I also would prefer to have all the files in one FC database. (14 files, all around 300 MB, total features 2648770).

I attempted processing one GML file on its own. I also realized that the GML files are all displaying in ArcCatalog as feature datasets. Each of the TOP10 GML files appears to contain the same feature classes (33 in total per feature dataset, or GML file in this case). However, when I processed one of the GML files into a geodatabase, only two of the feature classes from the GML file are present in the GDB. It said the tool was successful, but a warning within the messages reads:

  xfmap- Can't determine axis order from coordinate system or settings. Using assumed axis order of 1,2.

I'm not sure how to follow up here; does anyone have specific experience with Dutch TOP10 NL files?


texts: Duplicate print composer or use template?

texts:
I want to duplicate a print composer in a project to use with a different atlas coverage layer but I am having problems generating the atlas.

I have tried to both duplicate the print composer in the composer manager and loading from a template but neither method would generate the atlas. I'm getting the error message:

no matching atlas features found

I have used the same atlas settings successfully in another print composer (using the same data sets) so I know it should work.

Obviously, I could set up the print composer from scratch again but I want to keep the exact layout. It seems like QGIS is having a problem 'refreshing' the atlas coverage layer. Has anyone else encountered this? or have suggestions as to how I can get it to recognise the new atlas layer?


texts:
system description: Windows7 32 bıt OS, Grass GIS 7.0.0 Gnuplot 5.0

Solution attempt 01:

On the basis of following thread, commented out said path check lines from i.spectral.py script:
https://lists.osgeo.org/pipermail/grass-user/2006-June/034856.html

Solution attempt 02:

Set GRASS_GNUPLOT=C:\Program Files\gnuplot\bin (installed location of gnuplot.exe) throgh
settings-> Grass Working Environment-> Change Settings (g.gisenv)
based upon recommendations from the following discussion:
http://osgeo-org.1560.x6.nabble.com/Error-in-windows-version-of-i-spectral-td5018720.html

However, still unable to run i.spectral. Following error message appears in command output. Any workaround?

Command:

  i.spectral -g --verbose
  raster=ATMS_B1@Target,ATMS_B2@Target,ATMS_B3@Target
  coordinates=202334.820717,2645169.44223 output=SPCTRL

Error Message:

Traceback (most recent call last):
  File "C:\Program Files\GRASS GIS
7.0.0/scripts/i.spectral.py", line 241, in <module>
    main()
  File "C:\Program Files\GRASS GIS
7.0.0/scripts/i.spectral.py", line 234, in main
    draw_gnuplot(what, xlabels, output, img_fmt,
coord_legend)
  File "C:\Program Files\GRASS GIS
7.0.0/scripts/i.spectral.py", line 144, in draw_gnuplot
    grass.call(['gnuplot', plotfile])
  File "C:\Program Files\GRASS GIS
7.0.0\etc\python\grass\script\core.py", line 72, in call
    return Popen(*args, **kwargs).wait()
  File "C:\Program Files\GRASS GIS
7.0.0\etc\python\grass\script\core.py", line 56, in __init__
    raise OSError
OSError


texts: Iterate through each county in a shapefile, select the county and make a second selection by location - modelbuilder

texts:
The goal of my project is to create a table with two fields. The first field showing a county name and the second field listing every county that touches the boarder of the county listed in the first field.

I am attempting to do this in ModelBuilder or Python. In ModelBuilder, I use the Iterate Feature Selection Iterator to select each county in the shapefile. This creates two outputs, Value and Selected Features. I then use the Selected Features output to make my Select by Location. Next, I use the model only tool Get Field Value. The hope is that this will get the field values of the NAME field in my attribute table for all of the selected features. Finally, I am using the Calculate Field tool to add the names of all the boundary counties.

Problems:

The problem with using the Selected Features output is that it isn't actually just selecting a county it is creating a temporary feature layer of just that county. So when attempting to make a select by location there is nothing else to select from.
The problem with using the Get Field Value tool is that is only selected the first field value in the table, rather than selecting all the values in the NAME field.

Does anyone know how I can fix these issues? Is this process even possible in ModelBuilder?

Python Script:


texts: QGIS 2.0.1.3 (Dufour) and 2..2.0.1 (Valmiera) 64 Bit versions: uninstall not completely possible, crashes at exit

texts:
haven't tried the stackexchange page before, but since I can't find information on this problem: I can't completely uninstall QGIS, I want to do a complete, clean reinstall, but whenever I start the installer after having uninstalled QGIS it shows:

  QGIS Dufour or QGIS Valmiera already installed

So it seems the uninstall doesn't work properly.

I'm working on a Fujitsu Esprimo with an Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz (8 CPUs), ~3.4GHz, 32768MB RAM, Windows 7 Enterprise 64-bit (6.1, Build 7601) Service Pack 1.

I've recently tried installing OSGEO4W in addition to QGIS. Can this be the problem? Another thing that kept happening with all my QGIS installation from 2.0 upwards (I have only the newest one installed currently, 2.2 Valmiera) is crashes at exit ("QGIS.bin has stopped working). I'm having no other problems with it though.

Could you give me instructions on how to completely remove the program and do a clean reinstall? Thanks in advance!


texts: Node JS and Tilestrata to serve SRTM Tiles

texts:
I emailed one of the authors and it turns out that tilestrata-disk wasn't meant to be used as a "provider" in their system. No wonder! Here's what Brian Reavis said about how the cache works:

  A cache expects tiles to come from somewhere else. It's job is to wrap a tile source (provider) with some sort of storage and retrieval mechanism. The disk cache has its own file and directory structure that it expects, so using it to serve your pre-cut tiles would probably be a problem.

So, I'm not going to do that. Now I will work on serving with a more complete stack that includes Mapnik for rendering via tilestrata-mapnik. He did rewrite part of tilestrata-disk so that it can be used as a "provider" should anyone care to try that.

I'm stuck with a server that runs without errors, but offers no output.  I don't know what's wrong, but my guesses are something to do with the raster's xml stylesheet (used by Mapnik to render the tiles from their GTiff source) or the source layer's projection (EPSG:3857 - spherical mercator).

This is what I have so far in my server.js file:

var tilestrata = require('tilestrata');
var disk = require('tilestrata-disk');
var mapnik = require('tilestrata-mapnik');
var strata = tilestrata.createServer();

//define layers
strata.layer('srtm')
    .route('tile.png')
        .use(disk.cache({dir: './tiles/srtm/'}))
        .use(mapnik({
            xml: './styles/srtm.xml',
            tileSize: 256,
            scale: 1
        }));
strata.listen(8080);

I've gotten rid of all the errors that showed up due to bad XML, but I can't get it the tiles to serve.  I when I navigate to http://localhost:8080/0/0/0/tile.png per Brian's suggestions, I just get a message that reads: "The image 'http://localhost:8080/0/0/0/tile.png' cannot be displayed because it contains errors."


texts:
I follow the guide pgRouting 2.0 for Windows quick guide to try using pgRouting to do some simple routing task. I successfully installed Postgres 9.5, PostGIS 2.2 and pgRouting 2.1, and used my own shapefile to test (I could not download test data described in the guide).

As described in the guide, I used the following SQL to make the line vector file routable:

alter table public.tc_2000 add column source integer;
alter table public.tc_2000 add column target integer;
select pgr_createTopology('public.tc_2000', 0.0001, 'geom', 'gid');

and the resulting table looked fine:

I then used the following SQL code (modified from the guide) trying to route:

SELECT seq, id1 AS node, id2 AS edge, cost, geom
FROM pgr_dijkstra(
'SELECT id, source, target, st_length(geom) as cost FROM public.tc_2000',
1, 2, false, false
) as di
JOIN public.tc_2000 pt
ON di.id2 = pt.id ;

then the error poped out:

ERROR:  An expected column was not found in the query
HINT:  Please veryfy the column names: id, source, target, cost
CONTEXT:  PL/pgSQL function pgr_dijkstra(text,bigint,bigint,boolean,boolean) line 6 at assignment
********** Error **********

ERROR: An expected column was not found in the query
SQL state: P0001
Hint: Please veryfy the column names: id, source, target, cost
Context: PL/pgSQL function pgr_dijkstra(text,bigint,bigint,boolean,boolean) line 6 at assignment

I am new to SQL, and I don't know how to modify my code to get results.


texts: How to load a memory layer into map canvas an zoom to it with pyQGIS?

texts:
I would like to load a vector memory layer into the map canvas with pyQGIS. Background is, I create a memory buffer layer around point data. While iterating the points I would like to visualize the buffer around each point in the map canvas and if possible automatically zoom to the memory layer. Is that possible?

Thanks


texts: How to get lat/long of all adresses inside a polygon in google map?

texts:
I have several polygon shapes inside a City in Germany. I need to generate Point data of houses inside those Polygons. It is very time consuming to add all the building Locations manually. I am wondering, if it is possible to extract the Point for all address somehow from Google maps as XML or kml Format.

For your Kind information, I do not have any Point data. I need to generate it somehow using Google map or Google earth from geocoded address. That is now the main challenge for me.


texts: Getting features for one polygon into another, Geometries do not line up using ArcGIS Desktop?

texts:
I have:
 - Census Blocks
 - Taxi Zones

The taxi zone shares a rough boundary with census blocks, but there are times where there is overlap due to a rough outline(In image, red line is Taxi Zone, Grey is census block) I know you could go in and manually 'snap' them to one another, but that is extremely time-consuming. I have been attempting to do this using Unions/Integrate/Elimination/Dissolve and other tools, but there always seems to be too many polygons because they get split up by both the Census Block and Taxi Layer.


texts: PyQgis getting attribute values from specific features without StopIteration exception

texts:
I have been trying to access the values of the attributes of certain features in a vector layer. I have tried to to as they suggest in How to read the attribute values using PyQGIS? and Getting feature attribute description with pyqgis?, but whenever i try to define the feature with next(), i get StopIteration exception. This stops my program as if it was an error. I have been working on it for quite some time, but I haven't been able to find a way around it.

Is there a way to get attribute values from specific features of the vector layer without getting stoped by "StopIteration"?

I am working in QGIS 2.18

EDIT, Code snippit:

for i in range(0, len(vaiable)):
    request = QgsFeatureRequest().setFilterFid(i)
    iterator = layer.getFeatures(QgsFeatureRequest().setFilterFid(i))
    feature = next(iterator)
    #feature = layer.getFeatures(request).next()

I have tried with both iterator and request, with but StopIterator is raised when trying to assign "feature" in both cases


texts: Clustering and combining points in PostGIS

texts:
I am trying to cluster points near each other and group them into one coordinate based on the centroid of that cluster based on a bounding box.  I'm able to cluster them and create a group number for the clusters, but I'm having trouble grouping them as one coordinate per cluster.
This is what I have so far:
SELECT ST_ClusterDBSCAN(geom, eps := 0.5, minpoints := 2) over () AS cid
FROM locations
WHERE geom
    @
    ST_MakeEnvelope (
        -114.496858, 33.149466,
        -112.006293, 39.082191, -- bounding box
        4326)
GROUP BY(geom)

However, this doesn't seem to group them based on the cluster id (cid).  I also need to return  the total of points within that cluster (I looked into ST_NumGeometries for this, but can't figure out how to integrate it into the query) and just the lat and lng.


texts: Raster calculator giving different results in ArcPy than in ArcGIS?

texts:
I have two rasters representing travelling speed: one for road (with speed values for different types of roads and "0" values where there is no roads) and one for land-cover (with lower values, but no "0").

I would like to combine them in a single travelling speed raster this way: Keeping road values where there is a road (so where the road raster is not 0) and, elsewhere, having the value of land-cover.

I do it all the time with the raster calculator in ArcGIS with the following command: Con(roads > 0 , roads, landcover)

I'd like to insert it in an ArcPy script so I did it this way:

raster = arcpy.gp.RasterCalculator_sa("""Con(roads>landCover, roads, landCover)""", tmpPath + "/friction")

When I do it in ArcGIS, it works perfectly well but when I do it in my ArcPy code, it only gives me the roads value as a result but ignore the landcover values...

I can't find why.

What could cause this expression to work in ArcGIS but not in my ArcPy code?


texts:
I've added an wms operational layer, I wonder if it's possible to add identify support, any samples will be helpful

{
            type: 'wms',
            url: 'http://...../geoserver/service_name/wms?service=WMS&version=1.1.0&request=GetMap',
            title: 'WMS Name',
            spatialReference: {
                wkid: 4326
            },
            options: {
                id: 'title',
                opacity: 1.0,
                visible: true,
                format: "png",
                spatialReference: {
                    wkid: 4326
                },
                visibleLayers: ['NAMESPACE:LAYER_NAME']
            },
            layerControlLayerInfos: {
                layerGroup: 'layer group name'
            }
        }


texts: Counting points B downstream of point A along polyline using GRASS?

texts:
How do I determine how many barriers are downstream of each sampling location on a river network using GRASS?

Resources

I have ~20,000 points, soon to be ~600,000 points, falling on or near a river network of Europe. The points are sampling locations and barrier locations. I also have nodes for the river network. Network and node layers (ECRIN for those familiar with European river networks) have attributes relating to segment_length, distance2mouthLine, is_source, is_mouth, basin_id and segment_id.

What I've tried

Thus far I've tried v.distance in GRASS to calculate the distance going upstream for each point along a river segment which I then subtracted from the segments distance2mouthLine to provide a distance2mouthPoint value for each point. Using this I counted all of the points with a distance2mouthPointless than each point with the same basin_id. Unfortunately this included barriers on different river branches rather than those along the path from the point to the river mouth.

I have investigated Linear Referencing in ArcGIS, however the distance along the network is the same as that which I calculated in GRASS and therefore has the same limitations. Network Analyst in ArcGIS also provides a similar metric and I'm unaware of how to automate it rather than manually select each start and end point.

Any idea if this can be achieved in GRASS by identifying something like each source to mouth path and using v.distance as above?


texts: data driven pages, point to ddp index polygon

texts: I am trying to figure out how to convert points(currently my ddp driver) into polygons that cover the extent of my data frame. The points feature class are often in close proximity to each-other, and I want to make tiles that overlap, so that I can use them to produce annotations for each tile and then query these annotations with tile ID's to prevent overlapping labels from the tiles.  I am not sure of a quick way to accomplish this.  I have thought about just creating book marks for each DDP and running  a script that converts those to tiles, but is seems like there is probably an easier way.


texts: How to get all band values with Identify tool of ArcGIS for Desktop?

texts:
This question is similar to How to always see all band values with Identify Features of QGIS?, but is for ArcGIS.

I have a multiband raster image, with more than three bands, and I would like to use the identify tool to get all the band values, not just the three bands that are displaid.

Of course, the workaround is to add each individual band or several 3-bands composites, but then I need to have those bands visible as well (which slows down the visualization) or to identify all layers (which slow down the identify). Furthermore, in this case, I have to click in the identify window, which takes me some time.

As a more generic feature, it would be great to have on the fly pixel based map algebra in the identify tool. But this is one step further.


texts: How to join a table based on postcodes in QGIS?

texts:
I have 2 datasets, a shapefile with postcode polygons and an excel table (.xls) with postcodes and building count.

How to represent the excel table in the map?

I know that in ArcGIS it's possible to add a table and then do a spatial join but no idea how to do that on QGIS. Or even just use the geocoding option

I downloaded the spreadsheet plugin and did a join, but I still wonder if there would be a better way of doing this


texts: Adding unique value to an attribute table starting at a known value

texts: I have an attribute table with a unique values under polygon ID.  Last number is 679.  I am adding new 100 new attributes to this table and want to give these rows unique identifiers starting at 680. I do not was to change my old unique identifiers and I do not want to manually enter the new values.  What is the code to complete that?  I know there is on I have done it before.


texts: Save layer file with different Definition Query with ArcPy

texts:
I would like to save a new layer file with different definition query, but receive TypeError: 'unicode' object is not callable.

import arcpy

lyrFile = arcpy.mapping.Layer("\path to  layer file\ layerFile.lyr")
field = 'Date_1'
values = "'1950'"
query=str(field) + "=" + str(values)
for lyr in arcpy.mapping.ListLayers(lyrFile, "Name of Layer"):
    lyr.definitionQuery(query)
    lyr.saveACopy("new layerfile.lyr")

Error I receive is:

lyr.definitionQuery(query)
TypeError: 'unicode' object is not callable

How can I adjust the code so could save layer files with different definition queries?


texts: Loading selected records in postgis raster table into QGIS using DB Manager

texts:
I've recently started using (out-of-db) raster data within postgis.

The raster data is an output from a metocean forecast model, with new data available every hour.

I'm uploading each new raster to a database, using something like this.

raster2pgsql -a -s 4326 -I -F -R path/to/tif/latest_output.tif postgis_table | psql db

As you can see, the -a flag is appending each new raster to the table, so I end up with a table where each record is a raster (sorry about the table):

+------+---------------+--------------------+
| rid  |     rast      |      filename      |
+------+---------------+--------------------+
| 1    | RASTER        | raster_0000_hrs    |
| 2    | RASTER        | raster_0100_hrs    |
| 3    | RASTER        | raster_0200_hrs    |
| 4    | RASTER        | raster_0300_hrs    |
+------+---------------+--------------------+

While this is really useful for running queries, I'm struggling to view the raster data in QGIS via the DB Manager.

Is it possible to "filter" raster tables in the DB Manager to load in only a select few rasters (e.g. raster_0000_hrs only), rather than adding the entire table, which invetibaly crashes QGIS? Are there any alternative plugins that I could use instead?


texts: Calculate Number Of Pixels Visible From Pixel

texts:
I am looking for a way in which to take an input DEM and create an output raster whose pixel's value is the count of pixels "visible" from that location.

I use ArcGIS with Spatial Analyst, but also Manifold, GlobalMapper or GDAL.

How would I go about this task?


texts: python postgis query in QGIS problem with encode

texts:
I have the latest QGIS 2.18.6 and i try to use some postgis query using python in QGS python console.

field characters from mytable is my non-English language and in postgis have encode 'utf-8' and any time where I want to work In QGIS never haven't problem with characters encode.

but any time to use python and postgis I have any time problem with encode characters.

here the code :

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from qgis.PyQt.QtXml import QDomDocument
from qgis.utils import iface
import psycopg2

constring = u"dbname='database' host='postgre' port='port' user='username' password='code'"
conn=psycopg2.connect(constring)
checkCur = conn.cursor()
sqlCheck = "select * from public.mytable;"
checkCur.execute(sqlCheck)
sqlCheckRess = checkCur.fetchall()

mylist=[]

for i in sqlCheckRess :
    mylist.append(i)
checkCur.close()
conn.close()
print mylist

that work nice but that I take in print mylist is un-humanity language

print like this :

 u'2247_\u03a5\u06c6\u03b9\u04c3\u03c4\u03ac\u05bc\u03b5\u03bd\u03bf \u03c5\u03c0\u03cc\u03b2\u03b1


texts: Exporting 60 PDF files from each MXD in ArcPy?

texts:
I need to export 60 PDFs files so I am thinking if there is a way around to export 60 without having to open each MXD file at a time ?

I am using ArcGIS 10.1 SP1  Build 3143 .


texts: QGIS 3 imports DATE type colum in Oracle Database as QDate, losing hh:mm:ss

texts:
I am at the moment working on a lift for an QGIS 2 plugin to work with QGIS 3.

The plugin specifies a QgsDataSourceUri and uses it in an QgsVectorLayer.
In the specific table, laying in an Oracle database, there is a column which is not named 'date' but a german name.
The type of that column is DATE type, therefore containing year, month, day, hour, minute and second.

In QGIS 2.10 this column was automatically converted to an QDateTime object by the constructor.
In QGIS 3.4.2 however this column is converted to QDate, losing the hours, minutes and seconds by that.

In the documentation I cannot find any explanation on how to specify the way a column is converted to prevent that data loss.

Am I missing an configuration option that achieves to specify the conversion of a column?

The loading procedure looks like this:

def fetch_layer(self):
    uri = QgsDataSourceURI()
    uri.setConnection(self.host, self.port, self.db_name, self.usr, self.pswd)
    uri.setDataSource(self.schema, self.table, "SHAPE", "username=\'{}\'".format(self.usr))
    layer = QgsVectorLayer(uri.uri(), self.layername, "oracle")
    print(f"This should not be QDate : {type(next(layer.getFeatures()).attribute(self.date_attr))}")

For the case that this is indeed a bug, I created also a bug report at QGIS Tickts


texts: Geometry of layer disappears when creating scratch/memory layer with script

texts:
I'm trying to write a Python script (QGis 3.4.9) to create a group of scratch-layers (polygon, line and point). When it is inside the function the layer is created as it should, but when the function is done the geometry (and anything else it seems) is gone. The layer is there but no geometry. I'm clearly missing something?

def scratch_group():
    isthere = False
    sgroupname = "Scratch_group"
    root = QgsProject.instance().layerTreeRoot()
    for child in root.children():
      if isinstance(child, QgsLayerTreeGroup):
        if child.name().upper() == sgroupname.upper():
            isthere = True
            print("- group: " + child.name() + " already there")

if not isthere:
    root = QgsProject.instance().layerTreeRoot()
    sgroup = root.addGroup(sgroupname)

    stiep = "point"
    slayerp= QgsVectorLayer(stiep, "Scratch_point",  "memory")
    sgroup.insertChildNode(0, QgsLayerTreeLayer(slayerp))
    slayerp.startEditing()

QMessageBox.information(iface.mainWindow(),"Check","Geometry",QMessageBox.Ok)
return()
#---main---
scratch_group()
QMessageBox.information(iface.mainWindow(),"Check","If Geometry is gone",QMessageBox.Ok)

EDIT:
In the meantime i stumbled on Adding layer to group in layers panel using PyQGIS?
You have to create a clone of the layer put it in place and remove the original.
That partially solved it:

def create_memlayer_ingroup(stiep, sname, sgroup):
    vlayer = QgsVectorLayer(stiep, sname,"memory")
    QgsProject.instance().addMapLayer(vlayer)
    root = QgsProject.instance().layerTreeRoot()
    ilayer = root.findLayer(vlayer.id())
    clone = ilayer.clone()
    sgroup.insertChildNode(0, clone)
    print(ilayer)
    root.removeChildNode(ilayer)

def create_sgroup(gname):
    rt = QgsProject.instance().layerTreeRoot()
    return(rt.addGroup(gname))

memgroup = create_sgroup("Scratch_group2")
create_memlayer_ingroup("Polygon", "Scratch_Polygon", memgroup)
create_memlayer_ingroup("Linestring", "Scratch_line", memgroup)
create_memlayer_ingroup("Point", "Scratch_point", memgroup)

This works when the group is not placed at the top of the tree. When it is at top (if no other layer is opened), the second and third memlayer is opened twice within the group. So not perfect yet...


texts: Are there plans for improving performance of vector files in QGIS?

texts:
Does anybody know if there are any improvements planned pertaining to the performance of vector files in QGIS? Right now, a 1 GB file with contour lines needs several minutes to open.

In addition to that I have the suggestion that vector-files shouldn't be refreshed in whole when panning or zooming. This costs a lot of nerves and could be avoided through a cache.

Thank you.


texts: GeoPandas spatial join generates error in QGIS 3.8 Python console

texts:
I want to execute a spatial join using GeoPandas in the Python console QGIS 3.8 like this:

dfsjoin = gpd.sjoin(point, polygon)

but it doesn't work. It always shows the following error:

Exception ignored in: <function Handle.__del__ at 0x0000026EC6680BF8>
Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python37\lib\site-packages\rtree\index.py", line 895, in __del__
    self.destroy()
  File "C:\OSGEO4~1\apps\Python37\lib\site-packages\rtree\index.py", line 884, in destroy
    self._destroy(self._ptr)
OSError: exception: access violation writing 0x0000000000000001

But, if I run it in OSGeo4W Shell it runs correctly as following:

What is the different between python console in QGIS app and in OSGeo4W Shell?
And how can I fix it?


texts: Spatial Adjustment of Many Points in Geometric Network

texts:
I have been receiving surveyed water valve data and editing our point feature class accordingly.  The current procedure is to use the attribute transfer mapping tool and click from source layer to target layer for each surveyed point, with transfer geometry on. Instead of doing this manually I want to move all of the selected valve points from a feature class , which are connected by a geometric network, to the shapefile point which has been surveyed.

Is there any way to do this within an editing session?

I want to move the points with a geometric network connection to a new location. There is a unique field I can join the point feature class and shapefile to but I still don't know how to update the location for all the points at the same time.


texts: Create a Leaflet custom checkbox control

texts:
I want to create a custom checkbox control which will simply set a flag in jquery/javascript: if checked the flag = 'clustered' or if unchecked flag = 'unclustered'.
So far I have a control on the map but its not a checkbox and how do i get the state of the checkbox - if its checked/unchecked.

code:

function MapShowCommand() {
  alert("checked / unchecked"); //set flag
}

L.Control.Command = L.Control.extend({
    options: {
        position: 'topleft',
    },

    onAdd: function (map) {
        var controlDiv = L.DomUtil.create('div', 'leaflet-control-command');
        L.DomEvent
            .addListener(controlDiv, 'click', L.DomEvent.stopPropagation)
            .addListener(controlDiv, 'click', L.DomEvent.preventDefault)
        .addListener(controlDiv, 'click', function () { MapShowCommand(); });

        var controlUI = L.DomUtil.create('div', 'leaflet-control-command-interior', controlDiv);
        controlUI.title = 'Map Commands';
        return controlDiv;
            }
    });
var test = new L.Control.Command();
map.addControl(test);


texts: Merging lines that DON'T touch in PostGIS?

texts:
After a specific demand in my work that needed to create a single line from two untouched geometries (lines), I needed to develop a function and follow for contribution:

Passing as a parameter two separate non-touching geometries (lines), a single complete line is returned filling the empty space between the nearest ends.

Code:

CREATE OR REPLACE FUNCTION public.st_mergecloselines(
    geometry, geometry)
    RETURNS geometry
    LANGUAGE 'sql'
-- This function merge two lines that don't within (dont touching) returing a single multiline. By Emilson Ribeiro Neto - emilsonribeiro@hotmail.com
AS $BODY$

select ST_LineMerge(st_union(st_union($1,
    (case
            WHEN (st_distanceSphere(ST_StartPoint($1),ST_StartPoint($2)) < st_distanceSphere(ST_StartPoint($1),ST_EndPoint($2))
                and (st_distanceSphere(ST_StartPoint($1),ST_StartPoint($2)) < st_distanceSphere(ST_EndPoint($1),ST_StartPoint($2)))
                and (st_distanceSphere(ST_StartPoint($1),ST_StartPoint($2)) < st_distanceSphere(ST_EndPoint($1),ST_EndPoint($2)))
              ) THEN st_makeLine(ST_StartPoint($1),ST_StartPoint($2))

            WHEN   (st_distanceSphere(ST_EndPoint($1),ST_StartPoint($2)) < st_distanceSphere(ST_StartPoint($1),ST_EndPoint($2))
                and (st_distanceSphere(ST_EndPoint($1),ST_StartPoint($2)) < st_distanceSphere(ST_StartPoint($1),ST_StartPoint($2)))
                and (st_distanceSphere(ST_EndPoint($1),ST_StartPoint($2)) < st_distanceSphere(ST_EndPoint($1),ST_EndPoint($2)))
              ) THEN st_makeLine(ST_EndPoint($1),ST_StartPoint($2))

            WHEN   (st_distanceSphere(ST_StartPoint($1),ST_EndPoint($2)) < st_distanceSphere(ST_StartPoint($1),ST_StartPoint($2))
                and (st_distanceSphere(ST_StartPoint($1),ST_EndPoint($2)) < st_distanceSphere(ST_EndPoint($1),ST_StartPoint($2)))
                and (st_distanceSphere(ST_StartPoint($1),ST_EndPoint($2)) < st_distanceSphere(ST_EndPoint($1),ST_EndPoint($2)))
              ) THEN st_makeLine(ST_StartPoint($1),ST_EndPoint($2))
            else  st_makeLine(ST_EndPoint($1),ST_EndPoint($2))
            end)),$2))

$BODY$


texts: Issue converting DEM in ASCII with LO17 (apparent) projection

texts:
I'm a newbie where GIS is concerned, so please don't hesitate to point me to a manual.

We've got DEM files without any projection information. One is called 000426_LO17_200M_BD16_25.ORT and contains lines like the following:

-75925.00 2710900.00   834.88
-75925.00 2710800.00   832.98
-75925.00 2710700.00   830.77
-75925.00 2710600.00   832.38
-75925.00 2710500.00   830.83

So we think this is a LO17 projection, at a 200M DEM spacing, and we're guessing the 3rd column is elevation in m. Given that we know the data is from South Africa (or close to it), this seems about right.

We're trying to get this data into Google Fusion, but need to convert it into a format that Google Fusion will understand, or equivalently, into a gdal format that we can then convert to Google Fusion.

We've tried a gdal_grid conversion, but when attempting to import into Fusion, we get Fusion is complaining about no mapping for TMSO to proj.4

If anybody has any ideas or suggestions, it would be greatly appreciated.

All the best,
Craig


texts: Calculating bounding box by given center and scale in Android?

texts:
Given these conditions:

a scale like 1:50000
the center of the viewport is 100°E, 20°N
the size of the viewport is 400x600

How can I calculate the bounding box of the viewport?

The geographic coordinate system of the map is EPSG:4490.

We want to display them in a different projection like Mercator or latitude_longitude_projection(maybe this is the so-said un-projected).

The size of the viewport is in pixels.


texts: What is the most efficient method to iterate geoprocessing operations through an attribute table?

texts:
I want to calculate the clipped raster layer percentage for a shapefile consist of overlapping polygons, since ArcGIS doesn't support tabulate area for overlapping features, so I used a loop to process each of these overlapping feature one by one:

1: Create an independent shapefile for each polygon from the attribute table

2: Clip the raster layer using the single one polygon created from step 1.

3: Calculate the raster grid percentage by counting the numbers of raster grids in each categories.

I was able to code this up, and the script run well but slowly for a big dataset. I run this loop for 2000+ big overlapping polygons and it took hours. I tried to optimize this computation by:

1: Use arcpy.da.SearchCursor instead of arcpy.SearchCursor

2: Using in-memory workspace to storage the temperate files.

3: Use Parallel Processing: arcpy.env.parallelProcessingFactor = "100%"

These are all helpful but I still think can improve more. I am thinking about using python parallel processing tool such as map function. Could anyone suggest on how to apply the map function on this loop process?

The following is the skeleton of my script:

featureClass = "E:/ArcGIS/access_to_recreation/new.gdb/Portland_ME_service_2hr"
rows = arcpy.da.SearchCursor(featureClass, ['FacilityID'])

# Loop through attribute table

for row in rows:
    id = row[0]

 # Extract each row as a independent polygon layer throughout the entire table
    inputfc = "E:/ArcGIS/access_to_recreation/new.gdb/Portland_ME_service_2hr"
    outputfc = "in_memory\\feature"         # edit here!
    fieldname = "FacilityID"
    fieldvalue = id
    # construct the where clause by calling the helper function
    where_clause = buildWhereClause(inputfc, fieldname, fieldvalue)
    # execute Select analysis, create separate layers
    arcpy.Select_analysis(inputfc, outputfc, where_clause)

    # # Clip the raster data using each polygon
    inRaster = "E:/ArcGIS/PAD/For_Recreational.gdb/PAD_Parks_continental_final"
    inMaskData = "in_memory\\feature"  # edit here!

    # Execute ExtractByMask
    outExtractByMask = ExtractByMask(inRaster, inMaskData)
    outExtractByMask.save("in_memory\\extract_mask.img")

# build raster attribute table
    try:
        arcpy.BuildRasterAttributeTable_management("in_memory\\extract_mask.img", "Overwrite")

    # Count the 0, 1 pixels in the mask output
    raster_rows = arcpy.da.SearchCursor("in_memory\\extract_mask.img", ['OID', 'Value', 'Count'])

    try:
        for ra_row in raster_rows:
            if ra_row[0] == 0:
                blank = ra_row[2]
            else:
                value = ra_row[2]
        percentage = 100 * float(value)/(float(blank)+ float(value))
    except:
        percentage = 0

    Percentage.append(percentage)
    ID_list.append(id)

    # Delete intermediate files
    arcpy.Delete_management("in_memory")


texts: NAICS and SIC database by address/parcel number?

texts:
I have a large amount of addresses that I need to manually input each and every one of their building uses. There are well over a dozen of categories to choose from.

What kinds of databases should I turn to in order to expedite this process and is this publicly available information?


texts: I have followed this tutorial to create custom forms in QGIS (I work with version 2.18). However, I seem to be unable to control the shape and size of the form window. Thus, widgets are frequently out of scope or the window is oversized. Am I missing something ?


texts: Tried renaming and deleting the ESRI folder in C:\Users\Mus\AppData\Roaming, no luck... Using Windows 7 SP1 32 bit system. When I start ArcMap 10.3 it shows the the program is initializing and loading document... but nothing opens. But ArcCatalog is working fine. Any suggestion?


texts: Making selections from two tables based on the reference table in PostgreSQL/Postgis

texts:
Based on the address locations ad for the region of interest (as a reference), I am trying to make selections from two street layers street1 and street2. What I have tried is shown by the code below:

Select
    --from address table
    address_id,
    address_locations,
    --from street1
    foo,
    bar,
    --from street2
    alpha,
    beta
From
(
    Select
        ad.gid As address_id,
        ad.geom As address_locations,
        foo.street1 As foo,
        bar.street1 As bar,
        aplha.street2 As alpha,
        beta.street2 As beta,
        St_Distance(ad.geom, street1.geom) As d
    From
        public.ad, public.street1, public.street2
    Where
        ST_DWithin(ad.geom, street1.geom, 30.0)
        OR ST_DWithin(ad.geom, street1.geom, 50.0)
        AND ST_DWithin(ad.geom, street2.geom, 30.0)
        OR ST_DWithin(ad.geom, street2.geom, 50.0)
    Order By address_id, d
) As nested_query;

When I used only ad and street1, the execution time was 142 milli-seconds. When I also include street2, the query takes much longer. The query searches for streets within 30 meters and then within 50 meters from street1 and then does the same for street2 and returns desired values from the columns. This query is running (no syntax errors) but taking a very long time (Execution time: 30 minutes and counting...).

Can someone please suggest to improve the query and consequently the execution time? Is there any alternative way to do the same job? Address table has 87 rows while both street1 and street2 have 16,060 rows. All tables have spatial indexes. SRID of all tables has been set to be same i.e., 3044. Both streets tables are not exactly overlapping but more or less the same. They have been used because they have different columns regarding street data.


texts:
This is my first attempt to do a study on this topic.

I would like to get some advice on what data are required for doing flood mapping.

The area I have selected for this is a village which faces flooding every year.

I have got details about

Rainfall data
Soil type
Elevation of the area

And I have done some field surveying.

What other data are required?


texts: List index out of range when converting layers as a service definietion

texts:
In ArcGIS Pro python window:

I'm adding certain layers into a list. From that list, I am converting the list in a service definition draft. I use this same script all the time, sometimes I'll get this "List out of range" and sometimes I don't so I don't why I get this error. I get this error when it's trying to run the create web layer SSD function. I've tried to remove the "[0]" part (don't know what this means), but I get a different error because of it.

    import arcpy

arcpy.env.overwriteOutput = True

aprx = arcpy.mp.ArcGISProject('C:/naples/naples.aprx')
m = aprx.listMaps('Map')[0]

lyrs=[]
lyrs.append(m.listLayers('Installation')[0])
lyrs.append(m.listLayers('Building_Area')[0])
lyrs.append(m.listLayers('Storage Tank')[0])
lyrs.append(m.listLayers('Structure Area')[0])

arcpy.mp.CreateWebLayerSDDraft(lyrs, 'C:/naples/service_definition/naples_installation.sddraft', 'naples_installation', 'MY_HOSTED_SERVICES', 'FEATURE_ACCESS')


texts: OpenLayers polygon fill on IE has black shading

texts:
I have an odd problem on with OpenLayers on IE6, IE7 and IE8.  A newly drawn polygon has odd black shading.  In fact, even the line (stroke) has the same black outline.

Works great in Firefox, Chrome, Safari and even IE9.

I've tried everything from creating my own stylemap and symbolizers, to removing all styles and using default layer, and it still has a problem.

And to add to the puzzle, if you draw a second feature and intersect with the first, the intersected area has perfect colorization/rendering.

I've tried changing the basic style to have a different color (say blue) but then the color is again the base blue plus the black shading/fill attached to it.

Have run the basic drawing example, http://www.openlayers.org/dev/examples/draw-feature.html and it works just fine on all browsers.

Edit 1/10 5p PST

I'm removing my previous edit with "renderer" because it doesn't apply anymore.  however I was finally able to create a stripped-down sample that reproduce the error.

I wrote a page that accepts a URL param to enable the bug, and it'll be broken in IE8,7,and 6.  All you have to do is draw a polygon, and you'll see that it comes out black

http://dl.dropbox.com/u/47684/nc4/debug/RenderTest.html?break=1

And here's the source:

        <meta  http-equiv="X-UA-Compatible" content="IE=EmulateIE7" />

As you can see, just setting the layer OPACITY will break it. And this actually breaks on IE6 and 7. It works on 8, but I have compatibility mode set to 7, so yes that makes it break on 8.

Can someone from OL dev group comment why setting opacity is breaking?  I've traced the code down to Util.js and the modifyDOMElement function, and it looks to be ok.  But... ?


texts: Getting solid shape of dam lake in 3D by using QGIS?

texts:
I am working on a plugin. I will get a raster map of a dam lake for instance. When load that map I want to get a 3D model which fills the lake entirely. Then I will divide the lake on horizontal view and cross section (by giving segment number).

At the end of the day, I should get every box of this segmentation. I try to illustrate that above:

Is it possible that make in QGIS 3.x by using 3D attribute?


texts: Getting feature count of QgsVectorLayer using PyQGIS?

texts:
I'm using this code to get features from QgsVectorLayer

 QgsVectorDataProvider* provider=theVectorLayer->dataProvider();
            if(!provider)
            {
                return;
            }

            theVectorLayer->select(provider->attributeIndexes(),theVectorLayer->extent(),true,false);
    theVectorLayer->selectedFeatureCount();//here returns 0 but there is one polygon
    while(theVectorLayer->nextFeature(currentFeature))
            {
         ....//here Iget attributes and current geometry (one polygon)
            }

How do I get feature count?

Method featureCount() returns -1.


texts: I have a csv file which I converted to points.  It contains features that are a mix of lines and polygons.  Format is below.  Now I would like to use the point to line processing tool to generate the lines and polygons.  Not sure if my csv format needs to be changed but it does not separate each feature.


texts:
I would like to calculate a MSI index for the each image of the collection and I developed below presented script. But, somehow it is not working and it says "midir.divade is not a function". Can anyone please assist me in improving the below script?

var image = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
    .filterBounds(point)
    .filterDate('2015-08-01', '2015-12-31');

var addMSI = function(image) {
  var nir = image.select("B5");
  var midir = image.select("B6");
  var msi = midir.divade(nir);
  return image.addBands(msi);
};

var withMSI = image.map(addMSI);

print(withMSI);


texts: I've exported my timeline from google maps using their takeout function to perform some perfunctory analyses on them, mainly because i'm curious. I've exported in JSON and I'm struggling to find the format the time is measured in it looks like X (being seconds/ milliseconds/nanoseconds) since some date advancing forward, but I can't for the life of me find out what and from what date.


texts: How to apply convolutional neural network over a time series of Landsat images

texts:
I am implementing a fully convolutional neural network (specifically U-Net) over a time series of Landsat 8 images to predict land cover change. I am following the demo notebook for the predictions. However, this notebook only trains CNN over a single image (using median over multiple years) to predict single image output of impervious surface.

What I want is to input multiple months of imageries for training to predict the respective monthly land cover change (that varies spatio-temporally).

I have found Dataset.window() that can create windows from a time series of datasets (here), but have no idea how can it be implemented for my case of time series of Landsat imageries.

Any ideas or examples where this is implemented?


texts: Using Proj4js with OpenLayers to transform lat lng from 4326 to 27700 but coming out unchanged

texts:
Using Proj4js to transform lat lng from 4326 to 27700 but coming out unchanged. I am using OpenLayers 2.12 (OS API Leisure map) EPSG 27700, Proj4js (CDN)

I have added the definition to proj4

epsg4326 = new OpenLayers.Projection("EPSG:4326"); //WGS 1984 projection
        // To
        proj4.defs("EPSG:27700","+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs");
        EPSG27700 = new OpenLayers.Projection("EPSG:27700");

        var lonLat = new OpenLayers.LonLat(-0.1279688, 51.5077286).transform(epsg4326, EPSG27700);
        console.log(lonLat);

lonLat is coming out as  {lon: -0.1279688, lat: 51.5077286} so unchanged.

I think I'm missing a trick.


texts: How to configure ST_GEOMETRY Access for Oracle 64 bit?

texts:
I'm trying to configure the Oracle extproc to access the geodatabase with SQL, using the documentation here.

This is a Windows 2008 R2 server with 64 bit, and the Oracle server is also 64 bit.

I'm still getting the dreaded ORA-28595: Extproc agent : Invalid DLL Path when I run a spatial SQL query.

I tried to debug the issue by loading the dll using dependency walker, and saw that it was trying to load 64 bit dlls.

I know that SDE is still 32 bit, so I'm assuming that 32 bit-64 bit conflict is the cause of the problem.

I have configured this with a 32 bit Oracle RDBMS in the past.

Has anyone configured Oracle extproc to access the geodatabase with SQL with Oracle 64 bit?


texts: Selection by Location on Raster, Using a Vector Layer

texts:
I'm trying to identify raster cells by location to eventually convert them to vector polygons, and then add them to an existing vector layer. Below is a picture to illustrate the problem.

The green layer is raster and the purple is the polygon. Is it possible to add raster cells that intersect the polygon within a given distance (notice how there are portions of the raster off to the far left that would be problematic if I just converted to a polygon now and then did a union). Ultimately I want to identify areas for a given distance, then be able to convert to vector so I can union it with the existing purple polygon.

Is this possible? Everything I've searched for regarding selection of raster suggests using Zonal Stats - but that's not necessarily what I'm trying to accomplish here. Any suggestions or alternatives?


texts: WorldView 3 level OR2a processing - do I have data as DNs?

texts:
I have a World View 3 image which is level Ortho Ready Standard Imagery
(OR2A) processed and I wish to ascertain what this processing level means.

I need to work on surface reflectance and am aware that Digital Numbers (DNs) would need to be processed to Top of Atmosphere (ToA) radiance and then to surface reflectance.

However, being new to remote sensing the Digital Globe website information on what OR2a means has confused me (http://c-agg.org/cm_vault/files/docs/DigitalGlobe-Base-Product-FAQ.pdf) and I am left wondering - do I currently have data as DNs or as TOA reflectance?

Having checked, the pixel values for each band are integers ranging from 0 to 1227, so I do appear to have DNs!


texts: Exporting slope values from DEM-derived rasters to CSV from QGIS

texts:
I want to export the slope values with their corresponding latitude and longitude into a CSV file. There's a DEM of a region having an irregular shape (not square). I have to generate the slope rate map of the region using DEM then create grid points at 50m intervals and get the corresponding values of the slope at all the grid points with their latitudes and longitudes in a CSV format.
I had the slope raster map now but don't know how to proceed forward getting slope values in CSV format.


texts: How to create lines from points from different feature class ArcPy

texts:
I want to create lines based on two points:
1st line: point number 1 from pointFc, point number 1 from centFc
2nd line: point number 2 from pointFc, point number 2 from centFc
and so on

Where did I make mistake?

start_cursor = arcpy.SearchCursor(pointFc)
end_cursor = arcpy.SearchCursor(centFc)

desc = arcpy.Describe(pointFc)
desc2 = arcpy.Describe(centFc)
shapefieldname = desc.ShapeFieldName
OID_start = desc.OIDFieldName
OID_end = desc2.OIDFieldName

point = arcpy.Point()
array = arcpy.Array()
featureList = []

cursor = arcpy.InsertCursor(lineFc)
feat = cursor.newRow()

for start in start_cursor:
startFeature = start.getValue(shapefieldname)
start_id = start.getValue(OID_start)
pnt1 = startFeature.getPart()
for end in end_cursor:
    end_id = end.getValue(OID_end)
    endFeature = end.getValue(shapefieldname)
    pnt2 = endFeature.getPart()
    if start_id == end_id:
        point.X = pnt1.X
        point.Y = pnt1.Y
        array.add(point)
        point.X = pnt2.X
        point.Y = pnt2.Y
        array.add(point)
        polyline = arcpy.Polyline(array)
        array.removeAll()
        featureList.append(polyline)
        feat.shape = polyline
        cursor.insertRow(feat)

del feat, cursor, end_cursor, start_cursor


texts:
(This question was asked on another forum 3 months ago but received no response.  I'd really welcome some help getting over the problem)
My input map is a JPG scan of a not-to-scale map. I would like to warp it to EPSG:3857 so that it can forced to scale. The result is not warping as I'd expect.  Se images below.
The approach taken is to add 9 widely dispersed GCPs (ground control points) but no projection to a new TIF file. The GCPs have been added as X, Y, Longitide, Latitude, Z=0. gdal_edit does this step (from Python) without error and gdalinfo in.tif now shows:
    Size is 3421, 1925
Coordinate System is `'
GCP Projection =
GCP[  0]: Id=1, Info=
          (462,111) -> (-1.541,50.759,0)
GCP[  1]: Id=2, Info=
          (732,657) -> (-1.496,50.704,0)
GCP[  2]: Id=3, Info=
          (1998,1764) -> (-1.285,50.588,0)
GCP[  3]: Id=4, Info=
          (1911,687) -> (-1.291,50.701,0)
GCP[  4]: Id=5, Info=
          (1896,108) -> (-1.3,50.762,0)
GCP[  5]: Id=6, Info=
          (2691,360) -> (-1.158,50.732,0)
GCP[  6]: Id=7, Info=
          (2739,1098) -> (-1.149,50.659,0)
GCP[  7]: Id=8, Info=
          (3108,756) -> (-1.084,50.689,0)
GCP[  8]: Id=9, Info=
          (2163,1320) -> (-1.254,50.635,0)

The subsequent warp command is supposed to introduce a spatial reference system:
gdalwarp -of GTiff -t_srs EPSG:3857 in.tif out-warped-3857.tif

It runs, but the resulting map is stretched wrongly, from 3421x1925 to 3852x1381. Changing the output projection to EPSG:4326, the result is identical which I wouldn't expect.
There are three images to illustrate:
Reference map (OpenTopoMap)

Input map

Output map

Questions: Is this the best way to add projection and force a map to scale?
Why is the vertical dimension shrinking so much?


texts: Natural Neighbour to create DEM QGIS Hannover

texts:
I tried converting shapefile elevation coordinates using natural neighbour to create a DEM however I got the following error:

The following layers were not correctly generated.
C:/Users/stuar/AppData/Local/Temp/processing_eNJUbl/16ed8df2fe0e4bcfa26ec99fd0696461/TARGET_OUT_GRID.sdat
C:/Users/stuar/AppData/Local/Temp/processing_eNJUbl/5200f53077624943ac64aa1b4ef0a1c3/TARGET_OUT_GRID.sdat
You can check the 'Log Messages Panel' in QGIS main window to find
more information about the execution of the algorithm.

How can I fix this?


texts: Converting coordinates to lat lon in Madrid?

texts: I am using Python script for accessing this API for getting madrid traffic data. But I am really confused on how to get the geographic coordinates from this data. It has 2 fields st_x and st_y which I feel contain the information for the coordinates. I want to get those coordinates in lat long.


texts: Computing area of shared surface between adjacent 3D features in ArcGIS Pro?

texts:
ArcGIS Pro currently provides the Polygon Neighbors tool to compute the length of the shared border between features.
I want to do the same type of calculation for (closed) 3D features, i.e. I want to calculate the area of the shared surface between two adjacent 3D features. However, I can't find a tool in ArcGIS Pro for doing this.
Does this functionality exist in ArcGIS Pro?