texts: QGIS custom form and autocompleted data entering with "validation check"

texts:
In Qt-Designer I created a custom form where I can choose serveral species names in two different combo boxes. This combo boxes are linked via "value relation" to the accordant non spatial tables (gattung and art) in my SpatiaLite database. The whole species name consists of "gattung" and "art". Because these tables ("gattung" and "art") contain a lot of species names (latin) wich one can choose I think about an eased editing feature.

1) I'd like to generate an autocompletion when entering the first 3 characters into a line edit (instead of a combo box) for "gattung" so that on can choose only these species that beginns for example with "Que" from a list.

2) When I've chosen a name that matches with these 3 characters for example "Que" for "Quercus", in the second field "art" only valide names can be chosen. For example "petrea" ore "robur" as a subset of the whole "art-names" which one can choose. Valide names are "Quercus robur" or "Quercus petrea".
This could be done via line edit when entering 4 characters or with a combo box.

I thought about python and a lettercode (Quepetr for Quercus petrea) to solve this but doesn't know where to start and also where to enter the code for the widgets in QT-Designer.

I'm grateful for any help or idea!


texts:
I am trying to set up some GIS teaching with QGIS on a University network.

One issue is plugins, some I want to use are internal, like the Georeferencer, others are external contribs, like the Semi-Automated Classification plugin.

When installing a new plugin, does a user need admin rights? Does this vary depending on the type of plugin?  This needs to be run with students who have their own level of admin rights and it's hard to test before the actual first practical lab.


texts: Making Tiles - Varying the Zoom Level by Population Density

texts:
I want to create map tiles for all of the United States at high zoom level (12 or 13).  I am using TileMill and mapping census tracts (income).

However, I only need the high zoom level for populated areas.  How could I exclude the areas with low population?

I want to do this to minimize tile generation time and webhost disk space.

For instance, I'm going to treat Alaska separately as most of it has low population density.  I can also exclude oceans and large lakes.  I could probably exclude a large portion of WY, ND, SD, and several other western states.


texts: How to display points from Geoserver WFS on OpenLayers map

texts:
I'm try to display my points set from Geoserver in Openlayers map. I want to use WFS service, but at the moment I see one point, additionally in wrong place!
I want to display all at EPSG 4326. I transformed map to this coordinate system and my points on Geoserver (imported from PostGIS database) are at EPSG 4326.
I use Geoserver 2.15.1.

My code:

var wfs = new ol.layer.Vector({
    source: new ol.source.Vector({
        loader: function (extent) {
            $.ajax('http://localhost:8080/geoserver/pilgrimages/ows', {
                type: 'GET',
                data: {
                    service: 'WFS',
                    version: '1.1.0',
                    request: 'GetFeature',
                    typename: 'travels',
                    srsname: 'EPSG:4326',
                    outputFormat: 'application/json',
                    bbox: extent.join(',') + ',EPSG:4326'
                }
            }).done(function (response) {
                wfs
                .getSource()
                .addFeatures(new ol.format.GeoJSON()
                  .readFeatures(response));
            });
        },
        strategy: ol.loadingstrategy.bbox,
        projection: 'EPSG:4326'
    })
});

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
        source: new ol.source.OSM()
        }),
        wfs
    ],
    view: new ol.View({
        center: ol.proj.transform([0,0], 'EPSG:4326', 'EPSG:3857'),
        zoom: 6
    })
});

Can you help me?


texts: How are unsaved/temporary edits stored during versioned editing?

texts:
In ArcGIS Desktop, I can edit versioned, Oracle tables.
I generally edit the default version, and I've registered  the tables as versioned with the option to move edits to base.
I do things like:

Start an edit session
Insert, update and delete data
Undo and redo edits that were made in the edit session
Save the edits

The documentation tells me:

Edits are temporary until you choose to save and apply them
permanently to your data. You can also quit an edit session without
saving your changes.
When you save edits, you write them to the data source,
or a database.
Source: Starting an edit session section in  What is
editing?

Question:
The above documentation doesn't tell me how the unsaved/temporary edits are stored until they are saved. Are they stored in RAM, a scratch table, uncommitted rows, or somewhere else?


texts: QGIS Points in polygon - Difference between zero, and cells with no value?

texts:
I've carried out a points in polygon analysis. My output data contains the following results:

PolygonID    Point.Count
1            3
2            0
3
4            1
5            0
6            8
7
8            4

Evidently, the analysis is working in that it returned values for most polygons. However, some polygons have no values at all?

It is clear what polygons with a value of zero mean - there are no points within the polygon boundary in question.

However, why would some be returned with no values?

Can I treat these the same as zero? Or, is this perhaps a bug in QGIS?

EDIT:

Using QGIS Version 2.14

I took a CSV file with lat-long coordinates in WGS84 and then saved it as a shapefile.

This shapefile was then imported into a new project and re-saved as EPSG:27700.

It was then imported into a project with my polygons that were already in EPSG:27700, to carry out the points in polygon analysis.


texts: Measuring distance between certain relations in QGIS

texts:
I need to measure the direct distance between busstops in meters, not taking roads/networks into account.

I have the WGS 84 Coordinates + globalIDs for 28.000 pairs.

I already imported the Busstops from a csv (UTF-8) into a point layer.

My previous atempts:

Distance Matrix doesn´t work, because of the huge data size 28.000 x 28.000 is too much "waste" data and the operation fails...
"Point Connector" tool just creates "lines" without any measurments

The output should be 28.000 Lines with the distance in meters in a CSV, or just simply the measurment in m for every busstop relation.

Any ideas?


texts: How to add plugin icon to status bar using PyQGIS?

texts: I would like add my plugin icon to the status bar at the bottom of the screen. How would I do that? (QGIS API documentation on statusBar() seems to be leading to the QT docs and that site is down right now)


texts: Combining st_join with st_nn to include all points within, as well as within a given distance of a polygon

texts:
This is a follow up to my earlier question (Spatial join in R - Adding points to polygons with multiple corresponding points).

I have successfully joined a spatial points file to a polygon file in R using the st_join function within the sf package with more than one point being assigned to a polygon if necessary, duplicating rows but keeping all points which fall within a polygon.

st_join(polygons, points)

However I also need to join points which fall outside of the polygons but within 500m of a polygon to their nearest polygon. Points which are >500m away from a polygon can be discarded.

I thought that combining the above with st_nn from the nngeo package should work using the following:

st_join(polygons, points, join = st_nn, maxdist = 500)

However in this case only 1 point is assigned to a polygon, even if more than one point falls within a polygon or within 500m of a polygon. i.e. the rows are not duplicated.

Here is a screenshot of a sample of points and polygons:

And here is table showing how the points should be assigned to the polygons and how they have been assigned in the respective methods:

I find it a little strange that the second method does not keep the duplicates, even though it is based on the same function. Can anyone tell me what I'm doing wrong here?

Edit: I tried adjusting the k parameter but this simply joins the first points within the given distance up to the max number given and therefore can assign 1 point to 2 polygons. e.g.

st_join(polygons, points, join = st_nn, k = 10, maxdist = 500)

returns 5 points for polygon 89028 as there are 5 points within 500m, when in fact only 1 point should be returned (011-05-0529) as the other 4 points are already assigned to the polygon in which they fall. A point should only be assigned to one polygon.


texts:
I'm trying to identify crop row automatically (or as automated as possible). So far I can identify lines using supervised classification in QGIS (semiclassificationplugin) or Ecognition, I successfully extract the polygons but I can not turn the polygons into a single continuous line .
What I have:

What I want:

Today I saw this video https://www.youtube.com/watch?v=dwIDbXHS8nY. In the video, the software identifies the sugar cane and transforms into dots, then in a line connecting all the dots.

This software does what it's trying to do. I would like to know how this is done.


texts: ogr2gui geojson to csv do not export coordinates

texts:
I tried to convert ogr2gui to convert geojson files to csv. The geojson file contains point data only. The resulting csv only contains attribute data, no coordinates.

Is there a way to get the coordinates from geojson file using ogr2ogr?


texts: Where can I find equipment that will allow me to recreate a form of Google Streetview?

texts:
I am looking for more affordable equipment that would allow me to recreate the concept of Google's Streetview, as well as software that would make a comparable UI map (without it "belonging" to Google, ie. not having everything available via Google Earth). What I would like to create are scientific maps along rivers, streams, and wetlands, with sampling data (temp., PH, TDS, dissolved oxygen...) available to the user along the path in various increments. Essentially I am looking for camera equipment that can take 360° shots, that can record the images with X,Y values that can be input to a digital map and brought up on screen the way Street View looks. I know this is more of an digital imaging question, but it will be set up through a GIS system for user interface.
Thanks!


texts: I want to create a heatmap, I have a shapefile of us and geocoded 2000 addresses on it. Now I want to measure density of points and cluster by heatmap. I assume I need a longitude and latitude, but my addresses are based on street, city, state and country. How can I convert to x and y. Do I need to do that?


texts:
I am creating a soil property map. I am using SAGA to create my covariates.
For now, I am using SAGA wetness index for my 50 acre research area.
Couple of questions:

I was able to use the basic terrain tool and create my covariate (Sagawetness)
How to export it? I tried using TIFF and export raster but when I put in ArcMap, the coordinates are weird.
Also, when I load the sagawentess into R, it says the coordinates do not match.

I am using random forest to predict and need this covariate so I am lost.


texts: QGIS 3.0 migration - unable to install PyPDF2

texts:
I am migrating my Plugin to QGIS 3.0. Currently a bit stuck trying to install PyPDF2. I have been trying different things for the last couple of hours or more, without success.

I have got another version of Python installed (3.5, which I use for other things and 2.7...partly - I uninstalled it via Windows today although the folder is still there at C:\Python27. Probably I have to manually delete the remainder) Anyway I believe it should be quite possible to run more than one version. I have done it before.

QGIS uses Python 3.6, here: C:\Program Files\QGIS 3.0\apps\Python36

So in the standard scenario I would install PyPDF2 from OSGeo4W Shell. I run OSGeo4W Shell and type o-help which shows a list of available modules including pip3 (although not pip). I presume I should use pip3 to install PyPDF2. So I type this command:

pip3 install PyPDF2

I get the familiar message 'Fatal error in launcher: Unable to create process using '"''

I have trawled the forums for solutions to this but the ones I have tried have not worked for me. It sounds as if this should be fixed by upgrading pip, from what I understand(?)

If I run pip3 install --upgrade pip3 it doesn't work either. Same error, as might be expected

I have tried to run pip / pip3 from a command window opened from the path C:\Program Files\QGIS 3.0\apps\Python36\Scripts

If I do ./pip3 install PyPDF2 from a Powershell window there I get the same error again.

If I go up one level to the Python36 folder and run python -m pip3 install PyPDF2 I get the error No module named pip3

And so and so forth.

In an effort to make sure I am using the correct (QGIS) version of Python I added the '\QGIS 3.0\apps\Python36' folder to the PYTHON_HOME variable and added %PYTHON_HOME% to the Path variable.

How can I install PyPDF2??


texts: MySQL GeoSpatial function for ST_Distance_Spheroid? Metric used in return Type of GLength?

texts:
Is there a ST_Distance_Sphere, ST_Distance_Spheroid equivalent function in MySQL's GeoSpatial Support?
  If not, any way to emulate this?

  What is the metric used in return type of GLength (length of the LineString)?
  The manual for GLength says that GLength() is a nonstandard name. It corresponds to the OpenGIS Length() function. But I couldn't find any Docs for OpenGIS Length(). All I found was for ST_Length are they the same?

  Also for accuracy do I need to specify the SRID(4326)? How Do I do this in MySQL?


texts:
I have an excel (Fig) which has both X & Y in projected coordinates (WGS84/UTM). Is it possible to draw lines between the X & Y in QGIS ?

I've tried to create points based on X,From_y and X,To_y. I've created two points by just loading them using add delimiter layer (pasted To_Y below From_Y and loaded it). Then tried with Points 2 paths. But it didn't work out.

Anyone know how to fix this ?

All i need is traverse line between the points

Edit:

This is what i am looking for. I've coordinates(X1, a1 and X1, a2) for n number of points which has to be connected by a vertical line


texts: PostGIS spacial data type to represent geo-located polygons

texts:
I have a PostgreSQL database with circle shaped objects called 'spots'. Each of those objects has latitude, longitude and a variable radius, information that I need to define this circles. Given an specific coordinate, I select all the objects that contain that given point.

The thing is that I want to switch to PostGIS this stuff because of performance. I've started documenting about PostGIS but as a complete newbie I still have some doubts about data types.

What is the best way to store this kind of information in PostGIS? I don't need precision, because the radius can only have a limited lenght of 100km. I've thought about creating a polygon using ST_Buffer and save this into a geometry field in the 'spots' table. I also do not need an exact circle, a decagon is enough. The test code I'm using to know if this is meeting my needs is the following

create table spots (name varchar, geom geometry);

-- Insert a 2500m radius polygon with center in given coordenates
insert into spots
values (
    'Test',
    ST_Buffer(
        ST_SetSRID(ST_MakePoint(39.63467934415986, 2.6308822631835938), 4326),
        2500,
        4
    )
);

I have several doubts about this code. Once the geometric figure is created with ST_Buffer, is it still preserving it's position? If I provide the SRID code with given lat-lon coordinates, do I need to use ST_SetSRID or it's completely unnecessary?

The thing is that I've discovered that using ST_Contains is always returning true, with no matter about the lat-lon values I'm providing to the function.


texts: Getting "not data-defined" label properties

texts:
My plugin needs to read out the labeling of a feature layer. I can read via QgsPalLayerSettings the data-defined properties.

QgsDataDefined * dd = set->dataDefinedProperty(QgsPalLayerSettings::Size);
if (dd!=nullptr)
{
    if (dd->isActive())
    {
        if (dd->useExpression())
        {
            QString sExpr=set->dataDefinedProperty(QgsPalLayerSettings::Size)->expressionString(); }
        }
    }

But in most cases the user sets directly a value an then there is no dd! Surely it is very simple but I didn't find out.

How can I get label properties if they are not data-defined?

EDIT:

In between I found the member label (QgsLabel) of  QgsFeatureLayer, but the properties there are empty/standard.


texts:
I close a project with CARTO on the website (no through coding) but I don't see where and how to write a title on the export map.
It was possible with CartoDB but I didn't found solution with CARTO.
I didn't found any solution but looks crazy that none didn't think to write a "header title" to the map.


texts: Adjusting Points with Offset using ArcGIS Pro

texts:
I have a collection of GPS points that were collected along a road over the course of a few years.
Some years are right where they are suppose to be, but others are offset by 500 feet. Another year's data is offset by 700 feet.
Is there a way using ArcGIS Pro to adjust these points back onto the road without manually dragging them or snapping it directly along the centerline?
I would like to preserve the relative position of those points in relation to each other after they are moved.


texts: Lookup Projection Name from EPSG with pyproj

texts:
Is it possible to programmatically lookup a projection's name (e.g. Albers Equal-Area) from its EPSG using pyproj?

I see pyproj.pj_list, but it does not contain a mapping to the EPSG.


texts: How to specify transformation method when using ogr2ogr to reproject geometry?

texts:
I have a question regarding geometry transformation using OGR/GDAL (or perhaps another open-source solution).

Let's say I need to transform geometry from one coordinate system to another using specific transformation method.

For instance, I want to achieve this:

Input: EPSG:4269
Output: EPSG:4326
Transformation: EPSG:1188 - NAD83 to WGS 84 (1) - see https://epsg.io/1188

It seems to me, the ogr2ogr doesn't allow to specify the transformation mechanism.

Or am I missing out something?

Is there perhaps another tool I can use for this task?

This functionality is supported by ArcGIS tools, but I'm looking for an open-source solution.


texts: Adding labels to spatial grid in leaflet using R?

texts:
Problem:

I am trying to add grid cell labels to a leaflet map in R using shapefiles provided by TIGER/Line Census Data. The desired result is a leaflet map displaying the raster ID layer labelled from 1:ncell(r) that is clipped by the urban boundary of Albuquerque NM. I'm almost there, but I have been stuck trying to get the grid cell labels to show up when using leaflet.

The grid cells appear to be clickable, i.e. the cursor changes when hovering over the grid. However, no popup is displayed when clicking a cell.

The desired output does not have to be a popup. A static ID number centered within each grid cell would also be really helpful.

I am loosely following this tutorial.

Data Source:

2017 TIGER/Line® Shapefiles: Urban Areas

Code:

library(dplyr)
library(rgeos)
library(rgdal)
library(sp)
library(leaflet)
library(raster)
library(stringr)

# reading in shapefiles for the entire US
urban_areas <- readOGR(dsn = "./tl_2017_us_uac10.shp")

# pulling out only Albuquerque, NM shapefiles
albuq <- subset(urban_areas, str_detect(NAME10, "Albuquerque, NM"))

# determine the bounding box of the city boundary
bbox(albuq)
e <- extent(bbox(albuq))

# convert to a raster object
r <- raster(e)

# divide this raster into a 10 by 10 grid
dim(r) <- c(10, 10) ## spatially explicit cell sizes are ideal, but this works for now...
projection(r) <- crs(proj4string(albuq))

# add label ID to our grid cells
r <- setValues(r, 1:ncell(r))

# reconvert back to shapefile with the goal of creating a popup of the cell ID for each polygon
shape <- rasterToPolygons(r, dissolve = TRUE)

# clip the grid cells that contain the Albuquerque polygon
p <- shape[albuq, ]

# trim the grid perimeter to match the Albuquerque polygon
map <- gIntersection(p, albuq, byid = TRUE, drop_lower_td = TRUE)

# create a popup using the grid cell label layer
pop <- p@data$layer

# creating our map
leaflet(map) %>% addTiles %>%
  addPolygons(stroke = TRUE, opacity = 1, fillOpacity = 0.05, popup = pop
     smoothFactor = 0.5, color="black", fillColor = NULL, weight = 0.5)

Output


texts: Seeking Open Source implementation of Majority Filter?

texts:
I am performing postprocessing on Landsat GeoTIFF images and want to use the Majority Filter as mentioned in the book on Image Analysis, Classification and Change Detection in Remote Sensing and ArcGIS Pro documentation.

Is there a Majority filter implementation available which is open source?

Perhaps someone can suggest a filter in OpenCV which gives results similar to Majority Filtering.


texts: QGIS OSM Basemap with English place names?

texts: I am pretty new to QGiS and to cartographic software in particular. I wanted to create some maps of some cities in the Arab Gulf and wanted to use an Open Street Maps basemap. However, everything in the region is labelled in Arabic, which I can read but most of my target audience probably not. Is there a way to change the labels into English? Or do you know another good open source basemap for the region that I could/should use?


texts:
I am exploring QGIS 2.10.1-Pisa for the first time. I'm familiar with ArcMap, but I need to explore QGIS as an alternative.

I am having a hard time trying to figure out the units and scale situation. I've tried following a bunch of other forum posts with similar problems, but I'm still getting a wacky scale that doesn't seem to make sense.

I set the project CRS to State Plane Zone 3 ft (that's where my project and my data is), drag my aerial and shapefile in, but the scale is way wrong. The same data brought into ArcMap or AutoCAD looks right when I set the scale to 1:150. When I try the same in QGIS, the scale looks more like 1:20 or 1:40. I've tried different data with the same results.


texts: Available US State or National business (or employer) database lists?

texts:
This question dovetails with this post on geomarketing applications.

I am familiar with some of larger, paid-for business registry's like InfoUSA and Dunn & Bradstreet.  While the detailed information their business listings database include (like # of employees, annual sales, industry code, etc) is fantastic, if you are on a limited budget, sometimes the cost of purchasing the data outweighs the value you can get from it.  If you have geocoding capabilities in house, one could potentially create a reasonable surrogate by spatially-enabling (geocoding) a publicly available business listing on your own.

If you were looking to take on this type of project, where would you start your research? Do you know of US agencies at a State or National level that maintain this type of information, AND are publicly available for download or minimal purchase?


texts: Vendor WMS parameters not sent in GetMap requests via QGIS?

texts:
When adding WMS layers through Python plugin/script I run into issue with QGIS stripping of vendor parameters. I found following issue that was closed (https://issues.qgis.org/issues/12255). Proposed solution doesn't with checking settings in GUI "Ignore GetMap/GetTile URI reported in capabilities" is not good for my case, as I am doing it from Python script.

Is there way to activate this setting from Python script. Or is there any other work around?

Also the solution in this post (Non-standard WMS parameters not sent in GetMap requests via QGIS?) is not getting me anywhere, because WMS has to get & not encoded to parse it correctly.


texts: Empty output when clipping data in ArcGIS 10.3.1

texts:
I'm encountering a "empty output generated" warning when attempting to clip a data layer. The clip layer is a buffer drawn around roads while the input layer is a shapefile output from GeoMason. The input layer has been projected to be the same projection as the clip layer.

Executing: Clip input E:\LUP\scratch\input_Clip.shp #
Start Time: Sun May 28 23:08:15 2017
Reading Features...
Cracking Features...
WARNING 000117: Warning empty output generated.
Succeeded at Sun May 28 23:08:25 2017 (Elapsed Time: 9.79 seconds)

As you can see from the screenshot below, there are overlapping features between the input and the clip (red outline). I've also double checked the projection and extents and they match. What could be the cause of this error?


texts: How to open shapefile dataset using ArcGIS Pro SDK .NET?

texts:
I write AddIn for ArcGIS Pro using ArcGIS Pro sdk.

I have a shapefile on my disk(for example "D:\Data\shapefile.shp").

How can I open this featureclass?

For gdb geodatabase is very simple

using (var gdb = new Geodatabase(path))
{
 var dataset = gdb.OpenDataset();
}

but folder is not geodatabase and this code didn't work.


texts: The expression box has a handy "select from selection" button. Is there a way to emulate this with Python? I have the basic expressions down, but can't find a way to do this.


texts:
I have a toolbar (in ArcGIS) set up with two buttons using python addin manager, can someone help me make sure this works, when I run these buttons nothing happens.

import arcpy
import pythonaddins
import os

class ButtonClass1(object):
    """Implementation for ReconTB_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = True
    def onClick(self):
        os.startfile('subprocess.py')

class ButtonClass2(object):
    """Implementation for ReconTB_addin.button_1 (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = True
    def onClick(self):
        mxd = arcpy.mapping.MapDocument("Current")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        addLayer = arcpy.mapping.Layer(r"mypath/sample.shp")
        arcpy.mapping.AddLayer(df, addLayer, "TOP")


texts: Correcting slope of river bed on DEM using QGIS or ArcGIS

texts:
I have a DEM that has incorrect value of altitude/slope in the river bed area. This is a result of the interpolation of elevation value that are close to the river and creates lots of errors.

In fact (as you can see in the image above) in some points the bed river go up.

I would like to correct the value of the DTM in the river bed and set a plausible values.


texts: Selecting Shapefiles by number of rows using ModelBuilder?

texts:
I have several shapefiles in a folder and all of them are points. What I want is to choose those that have more than two records.

I am using the ArcGIS ModelBuilder. I made an iterate feature class and then I made a get count that give me the number of rows for each shapefile. From this step I cannot do anything else. But I think the solution is to use one tool that has a logical condition and so let me just obtain the shapefiles that has more than 2 rows.

Any suggestions?


texts: SAGA GIS raster calculator - QGIS scripting gt expression

texts:
I am trying to perform a simple raster calculator that takes input rasters and outputs a new raster with a constant value where values exceed a specified number in the input.

This is part of of a QGIS script and currently works using GDAL raster calculator. However I am having issues where python stops working when running the script using GDAL RC. Therefore i'd like to adjust my script to use the SAGA GIS equivalent to perform the same process:

processing.runalg('gdalogr:rastercalculator', raster,'1',None,'1',None,'1',None,'1',None,'1',None,'1','A>0',None,5,None,OutRaster)

My question is how can I change the line above in my script to perform the same calculation using the SAGA raster calculator? That is how can I write the expression which will output a new raster with a constant value where the input raster exceeds 0.


texts:
I'm quite new to QGIS and found the great plugin QChainage.

Simply put, it allows to create points along lines.
There are different modes available. Create points every X meters or Divide the line i X number of parts.

What I would need to do is to choose how many parts the line should be split into, according to a column named for example "Parts".

How could I achieve something like that (perhaps using the build in QGIS python console)?


texts: Importing CSV file without x or y lines but only with street addresses and not using MMQGIS?

texts:
Is there a way to import CSV file and display it on the map with only the street address and not x or y coordinates?

I don't want to use MMQGIS as it is limited.


texts:
I just noticed that in the arcpy module (Esri) there's a sys sub-module.  Does anyone know what the difference is in using one or the other?  I've been looking around and can't find anyone who uses arcpy.sys, so why would it be there?

edit: I wasn't clear, but arcpy.sys seems to be exactly the same as the python built-in sys module.


texts: Why are numeric values in attribute table not sorted correctly?

texts: I have a problem with the new QGIS version (2.16.0 or 2.16.1 Nodebo). I want to sort numeric values in the attribute table by clicking on the field name. But in contrast to older QGIs versions the attributes are not sorted correctly. It feels like QGIS recognizes it as text. but in the layer properties under fields it says the attributes are numeric (double or integer).


texts: How to get the radius of a "circle" in Openlayers

texts:
How do you get the radius of a circle in openlayers? I have searched quite a bit and cannot find any resource for this. please help.

Here is my code but I want the radius or circumference of the circle and not the area.

         // The Drag function is now made and called
         // Add Drag
     drag = new OpenLayers.Control.DragFeature(polygonLayer, {
     autoActivate: true,
     onComplete: displayArea
     });

    function displayArea(feature){
    var area = feature.geometry.getArea();

    alert(area);

    }


texts: error 'names' attribute [4] must be the same length as the vector [0] when extracting rasters with polygons

texts:
I have two rasters images and a base shapefile. I am using the extract function to extract pixel values that are positioned in the 809 polygons of my shapefiled. Both of the rasters are outside of my area of interest so they will return null values. They have the same CRS. However, one of them doesn't return me null values and gives me instead this error:

Error in names(x) <- value : 'names' attribute [4] must be the same length as the vector [0]

Code that results in null values:

shp=shapefile("E:/SUO-PUO/Ostalo/Sentinel/Potencijalna_rjesenja/R/Azo_plohe_33N_buffer30m")

image1=raster("G:/GIS/R20/S2A_USER_MSI_L2A_TL_SGS__20160317T153554_A003837_T32TQR_B05_20m.tif")

pixel_values1=(na.omit(extract(snimka1,azo_plohe, cellnumbers=TRUE, weights=TRUE)))

pixel_values1

:0 rows

Code that results in error:

  shp=shapefile("E:/SUO-PUO/Ostalo/Sentinel/Potencijalna_rjesenja/R/Azo_plohe_33N_buffer30m")

image2=raster("G:/GIS/R20/S2A_USER_MSI_L2A_TL_SGS__20160317T153554_A003837_T33TVJ_B04_20m.tif.tif")

pixel_values2=(na.omit(extract(snimka1,azo_plohe, cellnumbers=TRUE, weights=TRUE)))

pixel_values2

Error in names(x) <- value : 'names' attribute [4] must be the same length as the vector [0]
5.
`colnames<-`(`*tmp*`, value = nms)
4.
.local(x, y, ...)
3.
extract(snimka1, azo_plohe, cellnumbers = TRUE, weights = TRUE, df = TRUE)
2.
extract(snimka1, azo_plohe, cellnumbers = TRUE, weights = TRUE, df = TRUE)
1.
na.omit(extract(snimka1, azo_plohe, cellnumbers = TRUE, weights = TRUE, df = TRUE))

I do not see why the image gets the error. I do know that if I remove the

df=TRUE

I will not get the error and will get a list of 809 null values for each raster. But the thing is I really need to use (atleast I think I do) the df=TRUE which is crucial for the rest of my code (so if you could please take that as granted?).


texts: Converting each Rastercell to its own Point-Vector, Google Earth Engine

texts:
Related Post

My goal is to perfectly obtain one vector point per raster cell, located at the Center of the RasterCell
Running .reduceToVectors() with the scale argument retrieves vector points, but none fit the raster.

var coorsfeat = mean.reduceToVectors({
  reducer : null,
  geometry : geometry,
  scale:6800,
  //crs :mean.projection().getInfo().crs,
  //crsTransform :mean.projection().getInfo().transform,
  geometryType : 'centroid',
  labelProperty : 'Precip',
  eightConnected: false
})

Running with:

scale: mean.projection().nominalScale().getInfo() = 111319.49079327357 retrieves no features

Running with :

crs :mean.projection().getInfo().crs,
 crsTransform :mean.projection().getInfo().transform,

retrieves no Features

Is there a more elegant way of achieving the goal?

Link to GEE


texts:
Whenever i zoom in print composer it jumps either in or out in steps not smoothly. I can not zoom progressively neither can i choose an area by using the magnifier. I am working on a map of islands and they are either represented very small in a large mass of water or the outer islands are not visible anymore. It must be possible to exactly define the area that i want to represent in print composer.
How do i do this?


texts: QGIS 2.14.4 attribute table fields not visible in Statist plugin

texts:
I had one point shapefile to which I joined new numerical attributes by a csv text file.

I saved it by right clicking and everything works fine, I use it for calculations etc.

But when I use the "Statist" plugin, the new fields are not visible.
The fields are visible in the properties tab of the shapefile and they are noted as qlonglong type (Typename: Integer64, Lenght: 10, Precision: 0).

Any idea why this might happen?


texts: Returning unique value from table using ArcPy?

texts:
Using Arcpy in Python, I would like to return / print '0.15'  from the following table:

When I use the following code I get "(u'D', 0.20)" as my print statement in Python:

import arcpy
fc = "C:/Users/.../table_test2"
q = 'C'
with arcpy.da.SearchCursor(fc, ['Index_Score'], where_clause=q) as cursor:
    print row

(u'D', 0.20)

This is not what I am looking for as I am looking for the value '0.15'.

When I tried changing my query statement 'q' to:
    q = '"SA_ALL" = "C"'
I get the error:

" Traceback (most recent call last):
File "<pyshell#143>", line 2, in <module>
for row in cursor:
RuntimeError: An expected Field was not found or could not be retrieved properly.    [table_test2]"

Can someone help me and explain what I am doing wrong? Also what does the 'u'  mean in my first returned print statement?


texts: arcpy.PackageWorkspace doesn't keep reference to SDE enterprise geodatabase

texts:
I am publishing a geoprocessing service on ArcGIS Server v10.3.0 that updates data on an ArcSDE Enterprise Geodatabase.

I have the following at the beginning of my script:

arcpy.env.workspace = "Database Connections/TestDatabase@hq-gissql01.sde"

As discussed at Connecting Geoprocessing Service to ArcSDE?, when I publish the script as a geoprocessing service, that line gets changed to:

arcpy.env.workspace = arcpy.env.packageWorkspace

This is all well and good, but arcpy.env.packageWorkspace is simply the v101 folder for this geoprocessing task. It is not the actual .SDE connection file. I can confirm that the .SDE connection file does get moved to the v101 folder, but arcpy.env.packageWorkspace is simply a file path to the v101 folder--not that file.

Because of this, when I call arcpy.ListFeatureClasses() or arcpy.ListDatasets(), an empty list gets returned. From what I can tell this is because the workspace is set to a folder on the server's file system, and not an .SDE file.

How are .py files being re-written when they are published as geoprocessing services? What is the proper syntax to use when setting the workspace so that the published version of the script references an .SDE file and not simply a file folder?


texts: Efficient way of determining polygon intersection with bounding box

texts:
I have a set of polygons, and I need to check whether they are intersecting with a given bounding box(rectangle).  What I am doing is, I am taking every vertex of polygon and checking whether it's inside bounding box or not.

If yes
return true
else
Now I am taking every vertex(i.e 4 vertices) of my bounding box and checking  whether it is inside polygon or not,
using  the algorithm from http://assemblysys.com/php-point-in-polygon-algorithm/
if yes
return true
else
return false(box and polygon are not intersecting)

This way of approaching is taking too much time.  I want another algorithm which is faster than this.  I tried to search for an answer on Google but was not able to find anything.  I tried for finding code of mysql st_intersects() function on github but again I was unable to find that function code.

I know there are many algorithms but, because  I am new to this field I was unable to find algorithms,So I used the above approach.


texts: Python only using one core at a time to work?

texts: I noticed recently that python is only using one core at a time to work. Is there a way, maybe special builds or commands, to make it use all resources (second core)?


texts: How to make multipleQGIS mapbook pages in 1 pdf?

texts: Is it possible to export a mapbook in QGIS 2 with multiple pages all in one PDF?


texts: I have some polygons with attributes (in my case, some polygons have health establishments) and others without attributes (without health establishments). Is there a function that I can use in order to attribute health establishments to the polygons without it (using a distance approach)?


texts:
I have to merge almost 8000 rasters (TIF).

I tried using Mosaic To New Raster in ArcGIS 10.3, but it's processing really slowly and crashes quite often.

Is there a faster (or more stable) way to merge 8000 rasters?

I also have access to QGIS 2.10.1, but I'm willing to try any free software.


texts: Creating routable network table in PostGIS

texts:
I am trying to create the routable network table, as per http://underdark.wordpress.com/2011/02/07/a-beginners-guide-to-pgrouting/. Everytime I do it in PostGIS extension, QGIS crashes. The area is no so huge, consist of one British county. I'm also trying to do it with with pgAdminIII of PostgreSQL, it takes ages to do so. I do have the limited knowledge of queries, indexing etc. so I am not sure how to speed the process up. Is there something developed since 2011 (the last posts on that subject I found with this date) would enable me to do it efficiently? Thanks

UPDATE: (per @Pawel, Feb 19, 2013)

I have reinstalled PosrgreSQL and PostGIS. When I am trying to execute assign_vertex_id() I have: function assign_vertex_id(unknown, numeric, unknown, unknown) does not exist. What would be the problem?


texts: Creating a Predominance Visualization in QGIS

texts:
This is what I'd like to replicate. I'd like to show the dominant category for each region. Like the author, I also have 6 categories to show (languages).

Here is the article that is explaining his method. How do I go about replicating this for QGIS?


texts: How to measure Lake water surface temperature in GEE?

texts:
I need to calculate LSWT from MODIS in Google Earth Engine.
I tried to do my task with SST band of "NOAA AVHRR Pathfinder Version 5.3 Collated Global 4km Sea Surface Temperature" because my study area is Urmia Lake which is a large water body. but I think MODIS will be good.
Here is my code for NOAA:
https://code.earthengine.google.com/81a62ae32c3401f26ce2c176f3a385da


texts: QGIS modeler - use input filename in output filename

texts:
I need to use the input filename in the output filename in a QGIS model. I have spent a couple of hours trying to find a solution but can't.

If the input is contourXY.tab and I want output to be \refactored\contourXY_refactored.tab and \refactored\contour_5m\contourXY_refactored_5m.tab

This is easy in ArcGIS modeler (see here) but I need a QGIS process OR is this python solution the only option?


texts: OL3 bbox strategy on geoJSON vector layer?

texts:
I try to improve my current OL-3 webmap to load less payload using the bbox strategy. So I altered my code similar to the bbox examples out there:

overlayGroup.getLayers().extend([
    new ol.layer.Vector({
    title: 'Nodes',
    source : new ol.source.Vector({
      format: new ol.format.GeoJSON(),
      loader: function(extent, resolution, projection) {
          alert("lol:"+extend);
        var url = '/api/nodes?bbox=' + extent.join(',');
        $.ajax({
          url: url,
          success: function(data) {
            source.addFeatures(source.readFeatures(data));
          }
        });
      },
      projection: 'EPSG:3857',
      strategy: ol.loadingstrategy.bbox()
      }),
    style: createNodeStyle(),
    }),
    new ol.layer.Vector({
    ....

Unfortunatly nothing get's displayed, but also no warnings appear. I'cant't debug as no breakpoint triggers and Firebug doesn't list any transfer :-(
(If I remove the bbox specific stuff and witch back to a geoJSON source everything works fine)


texts: Why is ArcGIS Pro not symbolizing my entire range of values?

texts:
This may be an easy one but I am not a statistician so it's something I'd like an answer to.

I have a polygon feature class loaded in ArcGIS Pro that contains a field with values ranging from 1 to around 1200.  I am attempting to use graduated color symbology to display the relative difference in values based on that field.

However, no matter what way I symbolize it (Quantile is what I'm after, but have tried others), it will not pick up 1200 as the max value.  It lists 661 as the max value for the field, and then refuses to render values above that.

I have around 163,000 rows in the FC.  Of that, approximately 200 rows fall above 661.

How do I get ArcGIS Pro to make the maximum value for graduated colors, the actual maximum value of the field?

UPDATE: I tried a few more things:

Creating a completely new feature class and loading it from scratch
into ArcGIS Pro.  Didn't matter, still cut off at 661 when doing any
graduated symbology.
Loaded the same data into QGIS (3.4).  Upon setting it to graduated
symbology, it picked up the min/max values without issue and used
them properly for each of the preset classification schemes.

So it seems that Pro is doing some sort of automatic range boundary setting for large datasets; perhaps in an attempt to be "helpful."


texts:
I recently exported an ArcGIS model to python and was reviewing the script and am confused by the following line:

if TRU_W_DatableFeatures == '#' or not TRU_W_DatableFeatures:

What is the role of the '#'?


texts:
I have a raster layer that resulted from a series of calculations with

raster calculator (simple subtraction of two raster layers)
gdal_translate -ot Int16 -of GTiff ... (convert to Int16, as the result of above is in Float32)

however the gdal_translate command seems to forfeit the No-Data information. It ("correctly") assigns minmal values to the previously no-data cells (−32,768 in case of Int16), however doesn't store the No-Data information in the metadata (as seen in Properties of the raster layer, see picture).
Previous layer, result from raster calculator:

Problematic layer, result from gdal_translate with no No-Data value:

As a result, the values are displayed as regular raster cells by default. I'm aware that I can adjust styling and clip out of range values in styling options to make the "−32,768-No-Data"-cells invisible again. But I'd rather have the raster layer remember, which cell values correspond to "No-Data"!
I tried to assign it with the gdal command:
gdal_translate -a_nodata -32.768 -ot Int16 -of GTiff ...
however that shows no effect.
Can I manually set a No-Data-Value via Properties or some other tool?


texts: gdalwarp to clip a .tif: Works fine standalone, won't work in Python

texts:
I am trying to use the subprocess module to call gdal warp and clip a .tif using a shapefile. When I try the command in GDAL it works fine, but when I attempt to call it into a Python script it will not work, and merely crashes. Below is the code I'm using:

import os, gdal, subprocess, sys
cmd = ['gdalwarp','-cutline','C:/IDStudy/Input/Clip.shp','-crop_to_cutline','-dstalpha','C:/IDStudy/Input/Merged.tif','C:/IDstudy/Output/Clipped.tif']
proc = subprocess.Popen(cmd, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
stdout,stderr=proc.communicate()
exit_code=proc.wait()

if exit_code: #Oops, something went wrong!
    raise RuntimeError(stderr)
else: print stdout

This subprocess command works fine for merging/reprojecting but for some reason doesn't like to clip my raster...Any thoughts?


texts: Is it possible to add multiple labels for each element of a layer?

texts:
I want to add different labels for the elements of a layer.
Example:

One label with the number inside the marker.
Another label above the marker with the name of the element, different typography and size and with scale rule so that it is not displayed if the scale is above 20.000.
Etc..

As far as I have seen in my QGIS 2.6.1 I am only able to add one unique label.

Is there any way to add more than one?


texts:
New feature in QGIS 3.12 does not seem to work.

And it looks different from shown in changelog https://qgis.org/en/site/forusers/visualchangelog312/#feature-selective-masking

This is how it actually looks like

There is "Symbol layers to mask" part missing, with "This mask shape should be then selected as mask source in the Masks properties of a layer in order to enable masking" instead.

So I went to Mask properties, check boxes in Masked symbol layers and in Mask sources.

And nothing happened. The label looks exactly the same.
So, I oppened up Mask properties again, and found that all check-boxes were unchecked.

What am I doing wrong?


texts: Project loaded but Composer loading makes "Not Responding"

texts:
I am using QGIS 2.18.6 and have big project (200+layers) which doesn't have problem with loading, but there is Composer attached to it and it is loaded only to 50% and then it says "Not Responding". I never had this problem before. it actually started yesterday when I imported different Composer style - it crashed and from that moment I can't load my project.

Other Projects and attached Composers are working fine.

Is there any way how to remove composer from project before actually loading it?

I tried to Google where Composer files are saved but couldn't find it, even QGIS folders in my PC doesn't look like there is something saved.

I am using Windows 10 and have high spec PC so it can't be any issue with not enough hardware to load it.

Does anyone know how to fix it?


texts: First column of join attributes by location on QGIS

texts:
I am trying to join attributes by location using processing (I have many so I need to do them by batch)

processing.runalg('qgis:joinattributesbylocation', firstlayer, secondlayer,u'equals',0,0,'',1,layerOutput)

The first future are integers then the outputs became all negative numbers.
i.e. 4120034510 becomes -174932786

Could some one help me how to avoid this?  ( may be an easiest way to change first column to string , or better ideas)


texts: Leaflet JS add GeoJSON shape as a hole in a polygon

texts:
I have a dataset which already contains geojson. I can add it to my map with the following line of code.

L.geoJson(data.geojson).addTo(map);

And I get this

But I want to achieve this, where the surrounding area is covered up and my geojson shape is a hole in the middle of it.

According to the leaflet documentation

  You can also create a polygon with holes by passing an array of arrays of latlngs, with the first latlngs array representing the exterior ring while the remaining represent the holes inside.

So in theory something like:

var scotland = [[60,-13],[60,0],[50,4],[50,-13]];
L.Polygon([scotland,coordinatesOfShapeHole].addTo(map);

My problem is this function doesn't seem to like me passing in data.geojson or data.geojson.coordinates as the second parameter (coordinatesOfShapeHole);

Is there some way to convert data.geojson into a format the L.Polygon function is happy with?

Or alternatively do the same thing but with L.geoJson?


texts: pyinstaller3.4 package of pyqgis 3.7 standalone program failed import qgis.core due to missing PyQt5.QtPrintSupport

texts:
I have a qgis standalone program that is develop with python3.7/pyqt5.11 that come with the OSGEO4W's QGIS 3.4.4 (the latest LTR), it has the following imports:

import os, sys
import sip
from PyQt5 import QtCore, QtGui
from PyQt5.QtGui import *
from PyQt5.QtWidgets import QApplication
from qgis.core import QgsApplication, QgsVectorLayer, QgsProject
from qgis.gui import QgsProjectionSelectionDialog

Pyinstaller packaged it without problem, but when running its exe, it keeps saying:

Traceback (most recent call last):
  File "qgis_test.py", line 6, in <module>
    from qgis.core import QgsApplication, QgsVectorLayer, QgsProject
  File "C:\OSGEO4~1\apps\Python37\lib\site-packages\pyinstaller-3.4-py3.7.egg\Py
Installer\loader\pyimod03_importers.py", line 627, in exec_module
    exec(bytecode, module.__dict__)
  File "qgis\core\__init__.py", line 27, in <module>
ModuleNotFoundError: No module named 'PyQt5.QtPrintSupport'
[7724] Failed to execute script qgis_test

When I added the PyQt5.QtPrintSupport to hiddenimport, it will complain about not finding the next QtXXX library (e.g. QtSql, QtNetwork...) and on and on. This type of thing has been a reoccurring theme with each of the QGIS's new releases, really wish they can be consistent with these packaging conventions.

How do I deal with this?


texts: How to add all ".dat" files (raster) into QGIS, in a defined directory and defined legend (style), using python code

texts:
I'm using a QGIS 2.8.3 in Windows 10. I've added just one raster ("dat" file) to the QGIS, with a defined legend, using the following code. If I wanted to add all the "dat" files into a directory how could I correct my code to do that? Any idea?

The python script:

# coding: utf-8
#---------------------------------------------
from PyQt4 import QtGui
from PyQt4.QtCore import *
from PyQt4.QtXml import *
from PyQt4.QtCore import QSize
from PyQt4 import QtCore
from qgis import core, gui

lyr = QgsRasterLayer("C:/PATRICIA/IG/OPERACIONAL/FWI_"+file_date+".dat", "FWI_"+file_date)
s = QgsRasterShader()
c = QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
i = []

i.append(QgsColorRampShader.ColorRampItem(8.5, QtGui.QColor('#2b9e27'), '8.5'))
i.append(QgsColorRampShader.ColorRampItem(17.2, QtGui.QColor('#fee717'), '17.2'))
i.append(QgsColorRampShader.ColorRampItem(24.6, QtGui.QColor('#ec7000'), '24.6'))
i.append(QgsColorRampShader.ColorRampItem(38.3, QtGui.QColor('#cc292c'), '38.3'))
i.append(QgsColorRampShader.ColorRampItem(100, QtGui.QColor('#75263b'), '100'))

c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(lyr)


texts: Convert polygon to line using ArcGIS Desktop?

texts:
I have a polygon representing a wall.

How do I simplify this to a line representing the middle of the polygon using ArcGIS Desktop?


texts: I came across this Map Comparisons where they have used Waze map tiles. I searched around and couldn't find any official documentation for using Waze tiles in leaflet. I couldn't even find their privacy policy that concerns with usage of tiles. I'm retrieving tiles from http://livemap-tiles4.waze.com/tiles/{z}/{x}/{y}.png in leafletjs. Waze client apps are open source; so does it mean I can use waze tiles without concerning their policy?


texts: General Guidance on Integration of Sidewalks With Roads

texts:
I am trying to prepare a network dataset for network analysis. I have a road network that is pretty well constructed, and a pedestrian network that is not as well constructed.

The red Lines represent sidewalks and pedestrian paths, the blue lines represent roads.

Is there a way I can quickly and efficiently integrate the two to proper network connectivity standards? or will I have to go through by hand and connect everything manually? the sidewalks are often disconnected from each other where they intersect with roads, and the crosswalks are not modeled.

If there is no way to integrate them without a great deal of effort, would it be simpler just to indicate the presence of sidewalks as an attribute in the roads, and just connect the standalone paths? connecting every sidewalk and every crosswalk would be an enormous undertaking.


texts: Change CRS of layer, but keep location on map?

texts:
I have some layers I've inherited. They are all polygons. Most of them are in the correct CRS, which is EPSG 4326. Some of them are not, they are in EPSG 26918. They were natively created this way.

All the layers show on the map where they belong, I am assuming because the layers with the incorrect CRS are reprojected on the fly by QGIS.

I can change the CRS of the incorrect layers to EPSG 4326, but naturally, they now appear far away from where they should be on the map.

I feel like there must be an obvious answer I am missing here.
Is there an easy way to change the CRS of these layers from 26918 to 4326, and keep the location of the objects where they belong on the map?

A solution in QGIS (including pyqgis) or PostgreSQL would be great!


texts: Do Vector-tile file-names need zoom levels and if so, why?

texts:
I'm new to mapping and I'm looking to self-host some PBF-based vector tile files, mainly for a hybrid mobile app (Phonegap, HTML5, Openlayers or Mapbox-GL-JS).
If I need to generate or pre-render Vector-based PBF files (from OSM data), do I need to specify the 'XYZ'-Zoom level in the PBF filename (e.g. ''), If so, could anybody explain why?
To my uninitiated, noob instinct, the very nature of vector map files suggest that zooming/scaling shouldn't require multiple files for a given tile, I can understand that where bitmap/raster graphics are concerned.


texts:
I have a PostgreSQL/PostGIS table which stores data about all the warning signs around the world. The map data is provided to me by our vendor, and the data is such that no possible combination of columns from the table can be used to define a suitable primary key for the table.

My objective is to create a primary key and creating a new AUTO_INCREMENT-like column with plain integers is not going to cut it.

Since the only column with actually unique data is the geometry column provided to me, is there any way to use this data in order to develop/create a suitable primary key?


texts: QGIS plpgsql RAISE MESSAGEs do not appear in QGIS Log Messages Panel | PostGIS Tab

texts:
I refer to this post Show "Raise Notice" and "Raise Warning" in QGIS and to the herein referred issue https://issues.qgis.org/issues/12184#change-59414.

Using QGIS 2.18.18 LTR I try to display plpgsql RAISE MESSAGES from trigger functions and I think they should appear in the QGIS Log Messages Panel, PostGIS-Tab but they don't.

Is this the place for me to expect my plpgsql messages or do I miss something?


texts: How to add Unique ID number to a list if rows contain a consecutive sequence (Arcpy)

texts:
I am trying to use an Arcpy Python code to identify the Unique ID number for streets that contain consecutive addresses and ignore rows that do not. For example:

FULLST  |  LOR  |  HIR  |  UID
Gate    |   1   |  200  |   1
Gate    |  201  |  300  |   2
Smith   |   1   |  100  |   3
Smith   |  350  |  400  |   4
Smith   |  401  |  500  |   5

In this example I need to add UID numbers 1, 2, 4, and 5 to a list. The code I have so far (from a comment by @dslamb here) will grab 2 and 5, but not 1 and 4. For some reason it will not add the first street UID in the consecutive sequence.

streetlist = []
sc = arcpy.SearchCursor(STMS, None, None, "FULLST; LOR; HIR; MCN; OEN; UID", None)
previousSt = None
previousHIR = None
previousMCN = None
previousOEN = None
for row in sc:
    if not previousSt:
        previousSt = row.getValue("FULLST")
        previousHIR = row.getValue("HIR")
        previousMCN = row.getValue("MCN")
        previousOEN = row.getValue("OEN")
    else:
        if previousSt == row.getValue("FULLST") and previousMCN == row.getValue("MCN") and previousOEN == row.getValue("OEN"):
            if (row.getValue("LOR")-1) == previousHIR or row.getValue("LOR") < previousHIR:
                streetlist.append(row.getValue("UID"))
                previousHIR = row.getValue("HIR")
            else:
                previousHIR = row.getValue("HIR")
        else:
            previousSt = row.getValue("FULLST")
            previousHIR = row.getValue("HIR")
            previousMCN = row.getValue("MCN")
            previousOEN = row.getValue("OEN")
del previousHIR
del previousMCN
del previousOEN
del previousSt
del sc
del row

I am using ArcGIS 10.2, ArcInfo License level


texts:
In QGIS, I want to merge features lines that same attribute value and are adjacent. What characterizes a line to be adjacent to another is the final vertex of a line to be coincident with the initial vertex of another line.

The Dissolve command does not consider the adjacency criterion. The tool merged features that have equal attribute value. The figure shows an example, where numbers represent the IDs of each line. There are 3 lines that have the attribute value equal to 17. That way, if they have the same attribute and are adjacent then should be merged. There are also 2 lines that have the attribute value equal to 22, but they are not adjacent and so should not be merged.

The problem of the Dissolve tool in QGIS is that the criterion is based only on the value of the attribute. In this way, when executing the tool, the lines 22 will be unified creating multipart.

Another tool is the Singleparts to multipart which also does not consider adjacency. Would there be another way to merge lines in QGIS?

I know that in ArcGIS the Dissolve tool considers adjacency if "Create multipart features" option is unchecked. I'm looking for something similar.


texts: Installing GDAL/ogr2ogr in Ubuntu 18.04 returning Unmet dependency erorr

texts:
I am trying to install GDAL in Ubuntu 18.04 with the following commands
sudo add-apt-repository -y ppa:ubuntugis/ppa
sudo apt update
sudo apt install gdal-bin python-gdal python3-gdal

But it returning the following error
Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:

The following packages have unmet dependencies:
 gdal-bin : Depends: gdal-abi-2-4-0
            Depends: libgdal20 (>= 2.4.0) but it is not going to be installed
 python-gdal : Depends: python-numpy (>= 1:1.13.1) but it is not installable
               Depends: python-numpy-abi9 but it is not installable
               Depends: gdal-abi-2-4-0
               Depends: libgdal20 (>= 2.4.0) but it is not going to be installed
 python3-gdal : Depends: python3-numpy (>= 1:1.13.1) but it is not installable
                Depends: python3-numpy-abi9 but it is not installable
                Depends: gdal-abi-2-4-0
                Depends: libgdal20 (>= 2.4.0) but it is not going to be installed
E: Unable to correct problems, you have held broken packages.


texts: ERDAS Imagine: negative z-coordinates for automatically set tie points

texts:
I'm a new ERDAS Imagine user, so it's probably a beginner problem I'm facing: I get negative elevations in my DTM after automatic tie points setting..

I am geo-referencing a stereo-pair of ASTER images to a reference image (with x,y,z, coordinates) in order to generate a DTM from them. When I set the GCPs in the point measurement tool of LPS without the use of tie points, then triangulate, I can generate a quite useful DEM. However, if I use the option 'automatic tie generation' and check the coordinates after triangulation, I get, what seems, random elevations for the tie points (although x and y seem ok, both in the ASTER-images and the reference image): basically any value between -30'000m and +50'000m. Consequentially, the resulting DEM is useless (but this also confuses me somewhat, as I don't expect the 'DTM extraction' process to use the z-coordinates of the tie points, but only the ones of the GCPs).

Any clue what I am missing?


texts:
I am new to this stuff. Basically I am trying to show the data from my local database in the popup for some specific rivers, but I am having a hard time while enabling popup. I tried to find the solution but didn't find anything.

I am using Esri_Hydro_Reference_Overlay layer to show rivers in this map and I am getting this message:

the next message says

  This layer doesn't have any attribute data. You can, however, display attributes from a different feature layer in the pop-up.
  How do you want to specify what layer to use?

And it gives me two options

  Choose from my items

and

  Enter a URL for a layer

here is the Screenshot of it

Can anyone help me?


texts: Leaflet Basemap: Changing sources of basemap options in runtime?

texts:
I'm using Leaflet v 1.4.0. In LL, using L.control.layer, the basemap (and overlay) is assigned as follows:

$.fn.initBaseMaps = function (selectedIndexProjectionOption) {
    lyrNASAGIBSBlueMarble = L.tileLayer(nasagibs_static_templates[selectedIndexProjectionOption], nasagibs_static_template_options[selectedIndexProjectionOption]);

    lyrGebcoBathy = L.tileLayer.wms(gebcoWMSServers[selectedIndexProjectionOption], gebco_static_template_options[selectedIndexProjectionOption]);

    baseLayers = [lyrNASAGIBSBlueMarble, lyrGebcoBathy];

    // Define base map object
    objBasemaps = {
        "NASA Blue Marble": lyrNASAGIBSBlueMarble,
        "Gebco Bathymetry": lyrGebcoBathy
    };
    // Define overlays
    objOverlays = {};

    ctlLayers = L.control.layers(objBasemaps, objOverlays, {
        position: 'topright',
        unloadInvisibleTiles: true,
        collapsed: false // temporary
    }).addTo(map);

    //Initial choice
    map.addLayer(baseLayers[0]);
};
$.fn.initBaseMaps(selectedIndexProjectionOption);

This works fine. However, when I change certain 'display conditions', I need to redefine the source of the 2 base layers. How to achieve this?

Note that I am keeping the 2 radio buttons as placeholders and am trying to change the source of the basemaps in runtime. (Please note that I am not asking about how to toggle between radio choices - but change the underlying source of the basemaps).

It seems like even after programmatically removing the baselayers using map.removeLayer(layer), leaflet seems to be storing some internal reference to the old basemaps somewhere and forcing its display as well, along with the newer one. That is, I cannot completely get rid of the older basemaps.


texts: ESRI's new geoportal-server-catalog XSLT files

texts:
I am trying to implement the new geoportal 2.0.1 from ESRI's Github and I am having a lot of trouble finding the xslt stylesheets for many types, especially the arcgis-details.xslt
The file seems to be referenced to be in GEOPORTAL_INSTALLAITON/metadata/details/arcgis-details.xslt:

"arcgis": {
key: "arcgis",
evaluator: G.evaluators.arcgis,
interrogationXPath: "/metadata/Esri/ArcGISFormat",
identifier: "ArcGIS-Metadata",
detailsXslt: "metadata/details/arcgis-details.xslt",
}

But the files are not included in the repo and I can't find them anywhere.


texts: Speed up classification in field calculator

texts:
I have a shapefile with a little over one million features (but it could also be 10+ million).

I need to categorize them by an attribute they have. It is a float value, and the categories are:

1: attrib < 0.1
2: attrib >= 0.1 and < 0.3
3: attrib >= 0.3 and < 0.5
4: attrib >= 0.5

I want to add this category value in a new field.

Easy enough using the field calculator, I thought:

CASE
    WHEN attrib < 0.1 THEN 1
    WHEN attrib >= 0.1 AND attrib < 0.3 THEN 2
    WHEN attrib >= 0.3 AND attrib < 0.5 THEN 3
    WHEN attrib < 0.5 THEN 4
END

into a new field ('category', integer, length 1).

This process is very slow (about a minute per 1000 features).

If I just style the features according to the same rules (graduated, 4 classes), that happens in a manner of seconds and I get a fancy visualization, which I take to mean there has to be a faster way to do this.

So the question:
Why is the field calculator so much slower than the graduated style, and what can I do to speed this up?
I am not opposed or unfamiliar with pyqgis, so a python/processing script solution is just as fine. I attempted a small script that basically does the same as the field calculator above, but it has the same low performance:

##input_layer=vector
##General Tools=group
##Test=name

from qgis.core import *

layer = processing.getObject(input_layer)
layer.startEditing()
for i, f in enumerate(layer.getFeatures()):
    level = f['attrib']
    if level < 0.1:
        f['cat'] = 1
    elif level >= 0.1 and level < 0.3:
        f['cat'] = 2
    elif level >= 0.3 and level < 0.5:
        f['cat'] = 3
    elif level >= 0.5:
        f['cat'] = 4
    layer.updateFeature(f)
    if i % 1000 == 0: # Only added this to track how slow the script is
        print i
        progress.setText('Feature ' + str(i))
layer.commitChanges()

Running QGIS 2.18.15.


texts: How can I use shapefile or spatialite data for creating formatted PDF reports?

texts:
I am currently working on a landscape inventory in QGIS. I need to create a report based on the data we collected, as a pdf. The usual way in our group would be doing everything by hand. I figured out that using the atlas feature, I can create maps for each object, which I could use in a Word serial print. I can also export the dbf attribute table as a csv, and use it in the same word document as a data source, including paths to both photos of the objects and said maps. Additional information, like areas and lengths, would be a plus in some cases, but I can live without them. (I would probably populate the attribute table with this information for the reports.)

I wonder, however, how professionals would approach this. Any suggestions,as specific as possible? I  alsready pondered R and knitr, but I'm not sure if there is an easier way.


texts: Setting zoom level from table attributes using ArcMap

texts:
Is there anyway to set the zoom level after double-clicking on black arrow in table attributes?
It seems to default to 1:4,000. As I'm checking over hundreds of areas and it would be nice to be able to automatically zoom in (around 1:600) to each item  on the table without having to scroll the mouse wheel each time.


texts: Seeking resources for python scripting in QGIS?

texts:
Possible Duplicate:
  QGIS tutorials and resources

I would like to know regarding the available resources for starting to learn about python scripting for QGIS.


texts: Plotting more than two layers on cartodb map visualization

texts: I have a country table that I am using to plot 3 layers using the visualization wizard. I am using CATEGORY for the first layer and BUBBLE for the second. However when I add the third layer with the options available, it can not be distinctly represented. I am using a basic free license. I would like to know if the third layer can be represented distinctly say by using a patter for example.


texts:
I know that my question is similar to this one "QGIS Line feature reshape issue" but I'm using QGIS 3.0 and the issue is still actual.

When I reshape (add a missing part) of a line (1) it reshapes the other one as well (2) but I've selected to reshape only "left" one.

Does anyone know if it will be fixed in the near future?


texts: Converting MrSID files to other formats in ArcMap?

texts:
I have some large MRSID files of counties in TN. They are 3 banded, about 4 - 7 GB, and 3rd generation.

I've been in the process of converting them to 3 banded TIFF files using Export Raster Data.  I want to eventually discern agricultural use from them.  The conversion process normally takes over 2 hours.

Is there a better method than this?

I am using ArcGIS 10.1 for Desktop.


texts: Are there Python packages for Landsat 8 processing?

texts:
While there are a number of python projects for MODIS acquisition and processing (pyModis and get_modis), I have yet to encounter a toolset for automating Landsat 8 processing in Python. There's discussion of GRASS-based approaches, and of challenges related to automated downloading, but I'm curious if anyone's come across Python processing options. My functionality wish list would include items similar to what GRASS offers:

data acquisition and import
Auto-enhance colors, natural color composites
Calculate Top-of-Atmosphere Reflectance and band-6 Temperature
Haze removal
Atmospheric correction
Cloud identification
image classification
Export of results


texts: Seeking polygon data of wine (viticulture) areas and their main grape or varietal

texts:
Where is coordinate data on the geographic boundaries of wine (vitaculture) areas and their main type of wine (grape or varietal)? a GIS shapefile or geoJSON file format would be perfect.
Ideally I would like the polygon areas for wines of USA or Europe along with a few attributes corresponding to each area.
But even just data for at least one major country or varietal would be nice.
I can handle the geographic coordinates in any format.


texts: Searching on the fly within buffer in QGIS

texts:
I am using "Unique Values" widget with "Editable" option in QGIS 2.18.3, for my "Name" field in POIs layer to eliminate spellings mistake.

"Unique Values" widget fetches all the distinct values in a drop-down list which matches with your text string from the layer.

Now I need a solution to fetch values only within the specific buffer of a newly added point. I mean, when I add new feature point it should fetch values within 30 meters buffer of my point, matching with my text string.

I need this solution to eliminate duplication with existing data which is very common in our Data Entry case. Is it possible to use some "Expression" in "Constraint" option?


texts: Is there a list available of plugins compatible with QGIS 3.0 Girona?

texts:
On February 23rd, QGIS 3.0 Girona was officially released.
I read in the documentation that the release contains breaking changes that will affect some or all of the existing QGIS plugins.

I've been searching for a list that tells whether

a plugin is known as being compatible with 3.0
a plugin is known as not being compatible with 3.0 and planned to be made compatible by its developer
a plugin is known as not being compatible with 3.0 and abandoned by its developer
a plugin has an unknown status

Compatibility doesn't seem to be shown on https://plugins.qgis.org/plugins/.
Does anyone know of (an initiative to make) such a list?


texts: What is a good taxonomy or naming convention for files and folders containing GIS data?

texts:
My company has collected about 30 TBs of GIS data over the past 8 years, and I always find myself asking the following questions:

What type of data do we have for a given geographic area?
What are the details about that data (e.g., resolution in meters per pixel)?
Where does the data exist on the hard drive so I can actually use it?
Have we processed the data already, or is it in an unaltered form from the source?

Up to and including now, I have attempted to address these questions by devising an appropriate folder and file taxonomy/hierarchy. Does anyone have any ideas/suggestions on some understandable, perhaps even standard ways of organizing GIS data using files and folders?

I am also open to learning more about how using a database might benefit my company; we are software developers, not GIS experts, so I suspect we're quite a bit behind the curve on how best to approach the problem of storing/organizing GIS data for ease of use. I did see the question Best practices for managing geospatial data but was only able to draw marginal use out of the answers because I am so unfamiliar with geodatabases.

UPDATE: This last week I spent a fair bit of time reading about GIS databases, and got started familiarizing myself with PostGIS. Long-term, I think we will end up moving towards employment of a database plus metadata server as recommended by JasonBirch in Best practices for managing geospatial data.


texts: Modifying Score text in Locate Widget of Flexviewer 2.4?

texts:
I want to add text where it says "score" to make it more clear that this is an address matching score.

I am using the LocateWidget_US.xml and do not think I can do it here.  I think the place to do it may be in the LocateWidget.mxml in this line????

 locateResult.content = getDefaultString("resultScoreText") + " " +
 (addrCandidate.score % 1 == 0 ? addrCandidate.score : addrCandidate.score.toFixed(1));


texts: Raster data type as input in processing algorithms

texts:
Why the processing algorithms run just with input files located in QGIS desktop? Is there any way to work with directories?

I am developing a plugin and I am using processing algorithms but the number of input files doesn't give the possibility to work in QGIS desktop. Is there any way to work with directories?


texts: Service areas for points not directly on lines

texts:
I am using ArcMap 10.4 to find network distance service areas to train stations on a road or footpath (not rail!) network.

My problem is that Network Analyst locates all stations, but solving the analysis does not build service area polygons unless I manually move stations onto a road junction first. For example, the gap on the left here is for a station I didn't move:

From the initial road featurelayer, I ran the Integrate, Feature-to-line, and Split-line-at-point tools to get valid segments.

I created the Network Dataset to not model turns; use 'Any vertex' road connectivity; not model elevation; not establish driving directions; and build service area index.
On solving, I also tried increasing the search buffer in the service area properties as here: https://gis.stackexchange.com/a/171425/88199.

What crucial step am I missing?

I'm working through the problem on the desktop before moving on to a code solution. I'm also open to pgRouting or QGIS if anyone has suggestions using either of those.


texts: How to access an ArcGIS Server secure service from java

texts:
I want to call an ArcGIS Server service through Java. The service is locked and only the publisher and administrator can access it.
If I know the credentials, what sort of authentication should I use? (I have web adaptor installed)


texts:
I am looking for a robust way to determine if a geo point (long, lat) is located on a road
I am currently using openstreetmap overpass-api with the following query:
way[highway](around:3.0, {lat}, {lng});out;

And then I am looking for relevant values of "way" with "highway" from the following:
'motorway', 'trunk', 'primary', 'secondary', 'tertiary', 'unclassified', 'residential', 'service', 'tertiary_link', 'motorway_link', 'trunk_link', 'primary_link', 'secondary_link'
Yet I see that a 'road' is somehow not counted for its full width, but only its "middle spine", Here is an example of a grid with my classifier (Red means on a road, Green means not on a road):

I can estimate the road's width by its "lane" count, yet this is only an estimation, and not all roads has the "lane" tag.
What would be the best way to approach this problem?

My data is only relevant to the US

I can suffer some level of mistakes / errors

A batch / bulk solution is preferred
EDIT:

The picture demonstrates the problem, I can only be close to the road's spine, and I cannot tell if I am on a road if I am far from the spine,
It would be great if I could get the full-road-polygon and query over it


texts: How to fix QGIS crashes when using GRASS toolbox?

texts:
I've installed QGIS 2.0.1 on Ubuntu 12.04. exactly as explained on qgis.org . Unfortunately QGIS crashes when trying to use any GRASS model (r.aspect as example):
just opening QGIS, click on toolbox expand toolbox, click on GRASS function r.aspect (as example) and it hangs... no error message, no possibility to define data for input!

It's probably a local problem but any help is appreciated. I've had QGIS 1.7.4 installed prior.

Additionally I tried on a fresh ubuntu build (12.04) in a virtualbox. same results: when trying to us a function from Grass-toolbox it crashes QGIS :-(


texts: How to label adjacent pages in data driven pages

texts:
I'm trying to create a mapbook using ArcGIS 10.1 data driven pages. Its quite large and I want to create match lines or labels so that if you are looking at a page, there is text/labels showing you which page is adjacent to the top, bottom, left, right.

If I remember, correctly, this option was part of the data driven pages toolbar in 10.0, but it doesn't seem to be there now?


texts: Getting pyqgis geometry() into PostGIS database via psycopg2 - Python 3.7 QGIS 3

texts:
Using pyqgis in a QGIS 3 Plugin, I can get features from a layer and geometry from features.   Is there a way to pass the geometry into psycopg2 so that it will work in PostGIS?

I did find the PPyGIS http://www.fabianowski.eu/projects/ppygis/  but could not get it to work with Python 3.7

# Create table
    cursor.execute("CREATE TABLE structures (name VARCHAR(32) PRIMARY KEY, start VARCHAR(16), geom GEOMETRY);")

    self.get_layers()
    pipe_lyr = QgsProject.instance().mapLayersByName(self.pipe_layer)[0]
    if not pipe_lyr.isValid():
        self.popMsg('Layer is not valid')
        return
    pipes = pipe_lyr.getFeatures()
    for pipe in pipes:
        geom = pipe.geometry()
        # Insert some data into table
        cursor.execute("INSERT INTO structures (name, start, geom) VALUES (%s, %s, %s);", (str(pipe['name']), str(pipe['start']), geom))
        print( pipe['name'])

It errors out with

  TypeError: not all arguments converted during string formatting
  Traceback (most recent call last):
    File "C:/Users/Cary/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\invert_builder\invert_builder.py", line 288, in azure_db
      cursor.execute("INSERT INTO structures (name, start, geom) VALUES (%s, %s);", (str(pipe['name']), str(pipe['start']), geom))
  TypeError: not all arguments converted during string formatting<

This is a plugin for QGIS 3.


texts: Exiting Streetview and resetting tilt and compass in Google Earth

texts:
Is there any way to force no tilt and the compass to face north when exiting Streetview in Google Earth? I know I can just press r after leaving Streetview, but I wanted to know if there's a way to set some default behavior so that resetting isn't even necessary.
Possibly related question (with no answer) at Keeping Google Earth Pro always facing North Up?


texts: Cannot make model parameter from output destination in ModelBuilder?

texts:
I have a model that only takes in one user-defined object, all other inputs to the tool are (deliberately) hard-coded. In ModelBuilder, I wish to make two model parameters, the input object (which defines the extent of the final output raster) and the destination of the final output raster.

I am unable to successfully make a model parameter out of the output destination. I have set this as a model parameter but it does not appear as a parameter when I run the tool in ArcMap, and the option to clear the value in ModelBuilder is greyed out, see below screenshot:

What is the error in my approach?


texts: Building polygon from a single point (assumed to be in bottom right)

texts:
I have a set of points which represent the bottom left corner of a 50x50 meters grid. Could you please let me know if there is a way to build a polygon/square from this point to represent the 50x50 area?

There are over 10,000 points in the file.

My point data is in Latitude and Longitude.


texts: Drawing boundaries between discrete features based on attribute values in ArcGIS for Desktop?

texts:
There is a point feature class with a certain numeric attribute (let us say, populated places with population).

I wonder if there is, in ArcGIS, a way to draw polygons such that the sum of the attribute for all the points within each polygon does not exceed a certain value.

Does anyone know how this might be done?

I'm using ArcGIS 10.2 for Desktop.


texts:
I run myscript and get this error:

from qgis.core import * # @UnusedWildImport
from qgis.gui import * # @UnusedWildImport
from PyQt4.QtCore import * # @UnusedWildImport
from PyQt4.QtGui import * # @UnusedWildImport
# Import the code for the dialog
from delineationdialog import DelineationDialog
***(code continues)***

Error is:

Traceback (most recent call last):
  File "C:/Program Files/QGIS 2.18/apps/qgis/python/plugins/HidroTurk/HT_HJK/delineation.py", line 28, in <module>
    from delineationdialog import DelineationDialog
  File "C:\Program Files\QGIS 2.18\apps\qgis\python\plugins\HidroTurk\HT_HJK\delineationdialog.py", line 24, in <module>
    from ..UIs.ui_delineation import Ui_QgsOptionsBase
ValueError: Attempted relative import in non-package

here is the DelineationDialog:

from PyQt4 import QtGui
from ..UIs.ui_delineation import Ui_QgsOptionsBase

class DelineationDialog(QtGui.QDialog, Ui_QgsOptionsBase):
    """Set up dialog from designer."""
    def __init__(self):
        """Constructor."""
        QtGui.QDialog.__init__(self)
        # Set up the user interface from Designer.
        # After setupUI you can access any designer object by doing
        # self.<objectname>, and you can use autoconnect slots - see
        # http://qt-project.org/doc/qt-4.8/designer-using-a-ui-file.html
        # #widgets-and-dialogs-with-auto-connect
        self.setupUi(self)


texts: Possible to set color based on attribute while using 2.5D symbology?

texts:
I'm using 2.5D symbology to have my features display at different heights based on a certain attribute, but would also like the features displayed as different colors based on another attribute.

Is this possible without creating a 3d map view (in qgis3) or using the Qgis2Threejs plugin?

I'm using Qgis 3.4 but also have Qgis 2.18.


texts: Cursor missing in table view of ArcGIS for Desktop?

texts:
In QGIS attribute tables there is a blinking cursor when editing data in tables. This is a very helpful feature that makes the orientation esp. in longer attribute values quite easy, or even just possible as well:

In ArcGIS 10.4 I miss such cursor (not necessarily blinking...), and it is mostly impossible for me (apart from counting the cursor-left-keystrokes) to predict where a character is going to appear:

I must miss some settings (table appearance, font-size?).


texts: gdalwarp not clipping VRT to correct extent

texts:
I'm trying to clip a projected .vrt raster to an extent expressed in GCS coordinates and it's not doing what I would expect. Here is my command:

gdalwarp -te -133.08813624 56.69042536 -132.69405069 56.93121339 -te_srs EPSG:4326 -tr 30 30 input_DEM.vrt clipped_DEM.tif

However the output coordinates of the clipped_DEM.tif are:
-133.181430831 56.6899107291 -132.601585159 56.9311731495

This doesn't seem like much, but it's enough to completely throw off my process as I need a DEM to completely cover a certain extent (the polygon featured in the screenshot below). The clipped DEM, also seen below, should actually completely surround the polygon with a little bit (0.05 degrees) of padding. Instead it appears the clip's extent is stretched extra to the left and right, but cut off from the top and bottom. I've pulled the .vrt into Arc to look at it and it definitely has data (way) beyond the resulting extent.

I've used the same exact process before, only the raster I was clipping was a (projected) .tif, not a (projected) .vrt. As a result I'm thinking the latter has something to do with it. For now I can make my process work by adding 0.5 degrees to the edges instead of just 0.05 degrees, but I would like to know why this is occurring as I will need to further develop this process in the future.


texts:
I have got an existing .NET C# desktop application which already reads various GIS formats. I would now like to add support to read data from ArcGIS Online.

What options have I got to access information such as services, groups, search data etc from ArcGIS Online? Can I use ArcObjects, REST endpoints?


texts:
A Tilestache config file has a "png options" object that passes its config straight through to the Python Imaging Library (PIL). Setting optimize: true does some kind of optimization as the file is being written by PIL. I've found that using tools like optipng or pngout can yield further optimization.

Is there some hook I can use, or a way to tell Tilestache to run optipng and pngout on the tiles once they're cached?

I found a tool called picopt that does this sort of thing; I could run it on the cache folder, but I'd like to know how I can have Tilestache do this automatically.

Or, if I'm barking up the wrong tree, what's the best practice for optimizing map tiles?


texts: Analyze height data using a digital elevation map

texts:
i ran the following code on python 2.7:

import sys, struct
from osgeo import gdal
from osgeo import gdalconst

minLat = -48
maxLat = -33
minLong = 165
maxLong = 179

dataset = gdal.Open("C:\Python27\l10g.gz")
band = dataset.GetRasterBand(1)

t = dataset.GetGeoTransform()
success,tInverse = gdal.InvGeoTransform(t)
if not success:
    print("Failed!")
    sys.exit(1)

x1, y1= gdal.ApplyGeoTransform(tInverse, minLong, minLat)
x2, y2= gdal.ApplyGeoTransform(tInverse, maxLong, maxLat)

minX = int(min(x1, x2))
maxX = int(max(x1, x2))
minY = int(min(y1, y2))
maxY = int(max(y1, y2))

width = (maxX - minX) + 1
fmt = "<" + ("h"* width)

for y in range(minY, maxY+1):
    scanline = band.ReadRaster(minX, y, width, 1, width, 1,
                               gdalconst.GDT_Int16)
    values = struct.unpack(fmt, scanline)

    for values in values:
        try:
            histogram[value] += 1
        except KeyError:
            histogram[value] = 1
for height in sorted(histogram.keys()):
    print (height, hsitogram[height])

for which i am getting error:

Traceback (most recent call last):
  File "C:/Python27/demagain.py", line 34, in <module>
    values = struct.unpack(fmt, scanline)
error: unpack requires a string argument of length 3362

Why is this python unpack function giving an error?


texts: Create single-band thematic raster in Google Earth Engine from multi-band image

texts:
I am working on translating some models a colleague created in ERDAS Imagine into Google Earth Engine. The model specifics in Imagine are largely irrelevant but the general steps are:

Perform a tasseled cap transformation on a Landsat scene and retain brightness, greenness, and wetness.
Reclassify the 3-bands into bins (1 through 8) representing each pixel's deviation from the mean, scaled by the standard deviation of the whole image.
Combine these 3 bands into a single band image with 512 classes where the value represents a unique combination of brightness, greenness and wetness binned values. For instance, if in the 3-band image the brightness/greenness/wetness bin values of a pixel were 1/3/5, the single-band image should display a value of 135 for that pixel.

I have succeeded in getting the first two steps to work but now need to create the final image with the 512 discrete classes. How can I accomplish this in GEE? I have tried using .and(), .cat() and explored if I could convert the band values to strings prior to using the .add() function but have come up empty.
Code below is not very useful but will hopefully provide an example of the products I need to combine.
// binned, discrete brightness image
var brightness_bp8 = ee.Image('users/ryanpatrickrock/SAL_lcTimeSeries/raster/boulder_07132017_brightness_bp8')
                      .clip(ROI)
// Map.addLayer(brightness_bp8, {min: 1, max: 8}, 'brightness')

// binned, discrete greenness image
var greenness_bp8 = ee.Image('users/ryanpatrickrock/SAL_lcTimeSeries/raster/boulder_07132017_brightness_bp8')
                      .clip(ROI)
// Map.addLayer(greenness_bp8, {min: 1, max: 8}, 'greenness')

// binned, discrete wetness image
var wetness_bp8 = ee.Image('users/ryanpatrickrock/SAL_lcTimeSeries/raster/boulder_07132017_wetness_bp8')
                      .clip(ROI)
// Map.addLayer(wetness_bp8, {min: 1, max:8}, 'wetness')

// 3 band image of all binned values
var bp_cat = ee.Image.cat(brightness_bp8, greenness_bp8, wetness_bp8)
// print(bp_cat)
// Map.addLayer(bp_cat)


texts: Using Google Earth Engine to filter Landsat Daytime image?

texts:
How do I filter Landsat imageCollection by their time of collection?

For example, if they are daytime or nighttime, which might be specified as Ascending or Descending.

I don't know which metadata keyword to use for filtering out by that.


texts: Using "SELECT DISTINCT" in SQL Query for Labels?

texts:
Working with ArcGIS 10.2. I'm using GIS data downloaded from SHA and I'm using it to label roads. I have three different features I've selected to label (different kinds of roads). To keep it brief: I'm trying to prevent duplicate labels from popping up. If I use the SQL Query option I can choose which field to label by from the attributes table.

Some roads are labeled more than once, esp in terms of direction (duplicate labels pop up for roads that have an E and W direction etc). Via some quick googling I determined that a "SELECT DISTINCT" command exists for SQL but ArcGIS doesn't seem to recognize it? Would anyone know how to prevent duplicate labels using the options available to me? I accessed this pop up from Properties > Labels > SQL Query... for the GIS data.


texts:
Full-res, full-disk GOES-WEST images are available from
goes.gsfc.nasa.gov.  Is it possible to rig up a projection
string for these images?

We know the center pixel is at (0N, 135W) and the
satellite elevation is 35,786km.  ARC has a "View from Space"
variant of the orthographic projection.  Can this be modified
to provide Lat/Lon coordinates for pixels on the image?

It doesn't have to be perfect, but I'd like to overlay coastal
boundaries.

PS: Don't know offhand if 35Kkm is from sea-level or from the earth's center.


texts: Mixing externalGraphic and default point renderer on VectorLayer

texts: Is it possible to mix externalgraphics and the default point renderer for vector layers? I haven't been able to find an example that does this.


texts: How to increment a specific amount of features

texts:
I current have attribute bar all empty with null values.
I would like to fill that attribute starting from 10001, adding +1 to each of the features and having N before the number
So far I have come up with this.
'N'+to_string(1000+1)

This is only giving N10001 to all features.
How can I tell it that it needs to go alongside the column and increment it from the value the feature above has?
I have another attribute that could be used for comparison, attribute object id already has a number sequence for all those null features.


texts:
When programmatically creating / saving a QgsVectorLayer - implementation of the ExportSymbology functionality does not seem to exhibit the implied effects.

Does such a command not indicate that the layer style is restored when reopening or loading the shapefile at a later time (could be off on interpretation of this function)?

        writer = QgsVectorFileWriter(str(someFileName), "CP1250", some.fields, \
             QGis.WKBPoint, some.crs(), "ESRI Shapefile", \
             symbologyExport=QgsVectorFileWriter.FeatureSymbology)


texts: Automating procedure to add shapefile to multiple map templates using ArcPy?

texts:
I have a project that has me create a package of about a dozen maps for a client and I do this multiple times per week. I am trying to figure out if there is a way to automate this with python. I am hoping to create a script that takes a shapefile input from the user and inserts that shapefile into multiple mxd templates (FEMA, USGS, Soils, Aerials, etc).

I am very experienced with GIS but somewhat new to python.


texts: Creating Feature Class Based on Select By Attribute Using ArcPy?

texts:
Using ArcGIS 10.4.1 and Arcpy I am trying to select a shapefile by it attributes and create new feature class based on unique values. I am able to create new fs inside the Mapset.gdb/Linear but it is not producing based on "NEW_SELECTION". In Fact I am copying entire Roads.shp in each of feature classes.

How can I fix this?

import arcpy
from arcpy import env

def unique_values(table , field):
    with arcpy.da.SearchCursor(table, [field]) as cursor:
        return sorted({row[0] for row in cursor})

uniques = unique_values(r'D:/GIS/Roads.shp' , 'fclass')

for unique in uniques:
    arcpy.MakeFeatureLayer_management(r'D:/GIS/Roads.shp', unique)
    arcpy.SelectLayerByAttribute_management (unique, "NEW_SELECTION", " 'fclass' = '"+unique+"' ")
    arcpy.CopyFeatures_management(unique, r"D:/GIS/Mapset.gdb/Linear/"+unique+"")


texts:
I am wondering how web apps like Garmin connect or Strava calculate a "moving time" from a gps track.

For example I might record a GPS track on my Garmin GPS that has an elapsed time of 5:34:26 but they calculate the moving time as 2:33:39.

I can't seem to find an discussion or examples on how this is done, but maybe I'm searching for the wrong things.


texts: PyQGIS Set Extent for loaded layers - Refresh problem

texts:
I've written a plugin which loads a number of map layers (based on constraints) and I would then like to set the canvas extent to a selection of them.  Code very similar to that below works perfectly, except when the plugin imports layers into a project with no layers in. If you load a single layer before running the plugin then it performs fine.

Note that code similar to below is run after the plugin loads in a number of new layers using for example: iface.addVectorLayer(path, 'name', 'ogr').

I've tried mapCanvas.refresh() and .repaint() within the plugin but none seem to work - what am I missing?

extent = QgsRectangle()
extent.setMinimal()
layers = [ iface.mapCanvas().layers()[0], iface.mapCanvas().layers()[1] ]
for layer in layers:
    extent.combineExtentWith( layer.extent() )
extent.scale( 1.1 ) # Increase a bit the extent to make sure all geometries lie inside
iface.mapCanvas().setExtent( extent )
iface.mapCanvas().refresh()


texts: What's the best way to add a projection to GDAL?

texts:
I have a number of files in Maryland State-Plane (US foot) that have the following wkt:

PROJCS["NAD_1983_StatePlane_Maryland_FIPS_1900_Feet",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",38.3],
PARAMETER["standard_parallel_2",39.45],
PARAMETER["latitude_of_origin",37.66666666666666],
PARAMETER["central_meridian",-77],
PARAMETER["false_easting",1312333.333333333],
PARAMETER["false_northing",0],
UNIT["US survey foot",0.3048006096012192,AUTHORITY["EPSG","9003"]]]

OGRSpatialReference::autoIdentifyEPSG fails to recognize this.
As far as I can tell, the only difference in this wkt and the one for 2248 on spatialreference.org is some trailing significant digits, some omitted Authority codes, and a reversal in the order of the standard parallels.

[Edit] I've been looking at the pcs.csv file and while I don't know exactly how to read the columns, it appears that the parameters in the file for 2248 don't match the above wkt or the wkt on spatialreference.org. The standard parallels differ, as does the latitude of origin. This ups the ante on what the right action is.

What is the preferred method for adding this to the data files so that it will be recognized? Do I edit my pcs.csv file? The pcs.override.csv? Or something else? Changing the content of the file is impractical in this case.


texts:
I am using QGIS 3.4 and I have unintentionally labelled two lakes with the same name as they are listed as one feature in the attribute table. I cannot figure out how to change the feature attribute for one of them so that they are separated.

Any suggestions, or do I have to redraw these areas?


texts: Cursor help on large spatial query in PostgreSQL

texts:
I have two very large tables on which I want to perform a spatial query.

UPDATE table_a
    SET att_1 = b.att_1, att_2 = b.att_2
    FROM table_a AS a
        JOIN table_b AS b
            ON ST_Contains(a.geom,b.geom);

I am rather new to cursors, but here is my current attempt:

DECLARE curs1 CURSOR FOR
    SELECT att_1, att_2 FROM table_a;

OPEN curs1

FETCH NEXT FROM curs1
    WHILE @@FETCH_STATUS = 0
    BEGIN
        UPDATE table_a
        SET att_1 = b.att_1, att_2 = b.att_2
        FROM table_a AS a
            JOIN table_b AS b
                ON ST_Contains(a.geom,b.geom)
        WHERE CURRENT OF curs1

        FETCH NEXT FROM curs1
    END

CLOSE curs1
DEALLOCATE curs1 ;

My current error is:

ERROR: syntax error at or near "OPEN"

LINE 5: OPEN curs1

You may have surmised that all I want to do is take two attributes from my second table and add them to my first table where they meet spatially (table a is polygon data and table b is point data).

I am currently using Postgres 9.4 with PostGIS 1.2.

The two tables are massive, 25 million and 5 million respectively. This is why I need to go down the cursor route, the straight update query simply fails with such large datasets.


texts: I was plotting some points with the plot function in GeoPandas. I realized that when I plotted using Point.lat, Point.long, the points didn't fall at the right place but when I exchanged the coordinates to the format Point.long, Point.lat, the points fell at the right places. Can anyone explain this?


texts:
When I paste the following on EE online code editor then if offers me the option to convert to an import record.

var l8sr = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");
var roi = /* color: #d63000 */ee.Geometry.Point([-123.1231231231234, 38.1231231231234]);

One disadvantage I found on converting it to import record is that quick copy (Crtl+C) of code do not get those imports.

Instead, I can Ctrl+C them if I leave them as normal code.

And I do not actually know the advantage of having them as imported records. May I ask what's the advantage?