texts: How to fix the label position after using automated placement in QGIS

texts:
I'm making a city map using QGIS 3.14 and I want to label every street (that's around 5000 label, the map is an A0 size at 1/12 000 scale). There are no other label than the street name and all label came from the same layer. After setting rule based labeling and fine tuning all setting for all rule I got a pretty good placement for maybe 80% of the street name but there are still area where I need to manually adjust the positions to get a printable result.
The problem arise when a manually move a label (I use the labeling toolbar "move" and "rotate" label button with data defined layer fields, I do not use the auxiliary storage).
Each time I manually move a label, the labeling engine recalculate the position of neighboring labels and an area that start as all good but one label end as only one label good and all labels around to manually replace.
In a similar way when panning the map the label engine recalculate the label position based on the part of the feature that are visible so even slightly panning the map result in different a label placement.
As each time I zoom, pan or adjust one label position all other visible labels change their positions I can't evaluate if my labels placements are globally satisfying or not and I can't predict how the labels will be placed in the exported map.
To try to solve the problem I set the "Allow truncated labels on edges of map" and lock the scale of the map and use the magnifier to zoom on the map but that doesn't make any difference.
I also try to pin all labels, but pinning the labels reset the rotation to 0 and all labels became horizontal...
So my question is are there a way to fix the label position set by the automated labeling engine ? (maybe by using an expression to fill the label_x, label_y and label_rotation field with the value set by the labeling engine) or are there any other way to get good placement on very densely packed label without manually placing all label ?


texts:
I wish to display custom colours for the pixel values in my raster. How can I create a custom colormap (cmap)?
In the example below, I have specified "magma" as the colour map. However, I wish to define my own colors for specific range/values.
import rasterio
from rasterio.plot import show,show_hist

if __name__ == "__main__":

    input_Dir = 'generated_46.tif'
    src = rasterio.open(input_Dir)
    show(src,cmap="magma")

Here the result of the above code


texts:
In this map I am looking to display each image as a thumbnail from the GeoJson file.

For the moment I am using this code:

var myIcon = L.icon({
      iconUrl: '',
      iconSize: [40, 40]
      });

var geojsonLayerP = new L.GeoJSON.AJAX("https://raw.githubusercontent.com/marinespatialviz/marinespatialviz/master/litter.geojson", {
         pointToLayer: function(feature, yx) {
           var marker = L.marker(yx, {icon: myIcon});
           marker.bindTooltip('
' + feature.properties.description + '
')
           return marker;
         }
       }).addTo(map);

But this does't give me the image as an icon. I know what I am doing wrong here, and I am looking for the right method to display each image as a thumbnail at the right location from the GeoJson file.


texts: How to use SQL queries in arcpy.SelectLayerByAttribute?

texts:
I'm trying to select features from a layer that I know exist. I'm using SelectLayerByAttribute_management in arcpy.
I have a layer that is properly defined earlier in the code named "tourmaline". I'm trying to query this layer to return all features with value of "Occurrence" in the field STATUS__NEW.

This is what I have:

where = "Status__New = 'Occurrence'"
arcpy.SelectLayerByAttribute_management(tourmaline, "", where)

No features on the map are selected when this code is run.
I've tried variations of the where clause using double quotes and square brackets around field name, but none of them work.

Using ArcMap's Select Layer By Attributes tool to select the same features using the following SQL statement selects the features correctly.

SELECT * FROM tourmaline WHERE:
Status__New = 'Occurrence'

Why is this happening? Why is my SelectLayerByAttribute statement wrong?


texts: Add layer to mxd files ignoring specific item- arcpy

texts:
I try to add layer (abc) to several mxd files with arcpy and ignoring mxd file which called "land use" . I work with arcview 10.3 :

import arcpy,os,sys,string,datetime
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project"
Layer1 = arcpy.mapping.Layer(r"C:\Project\layers\abc.lyr")
for mxd in arcpy.ListFiles("*.mxd"):
    if mxd == 'land use':
        continue
        print 'land use no encountered'
    print mxd
    mapdoc = arcpy.mapping.MapDocument(r"C:\Project\\" + mxd)
    df = arcpy.mapping.ListDataFrames(mapdoc, "Layers")[0]
    arcpy.mapping.AddLayer(df ,Layer1, "TOP") # TOP \ BOTTOM \ AUTO_ARRANGE
    print 'AddLayer'
    mapdoc.save()
del mxd

as a result, arcpy add the layer to all mxd files and this is a bad result:

>>>
land use 2.mxd
AddLayer
land use 3.mxd
AddLayer
land use.mxd
AddLayer


texts: Checking if property exists using ArcPy?

texts:
I'm new in ArcGIS and ArcPy.

I've got f.e. a Layer called "Flurnamen_5000". How can I find out (arcpy) if there exists a "Text Renderers" in Properties --> Symbologie --> Darstellung (look at the printsceen)?

If there is no Text Renderers in the Layer. I'd like to remove the Layer.
I started with:

>>> mxd = arcpy.mapping.MapDocument("CURRENT")
>>> lyr = arcpy.mapping.ListLayers(mxd, "Flurnamen_5000")

but I have no idea, how to check if a "Text Renderers" exists.

FeatureType is simple.


texts: QuickOSM / QGIS exported polygon not matching OSM

texts:
I used QuickOSM / QGIS to grab some place data but the polygon / coordinates are not matching the original on OSM (nominatim). This is "Olive Garden" in Schaumburg IL. On OSM is is the below left image. After downloading / exporting using QuickOSM / QGIS it is the right image. I've not changed any setting. CRS is default 4326 W84.

On OSM the polygon is starting as: [-88.0332509,42.0502153]...
The exported QuickOSM / QGIS one is starting as: [-88.0332143, 42.050252]...

How to export this to match the original coordinate at OSM?


texts: Error when creating seach filter from geojson

texts:
I tried to read geojson object from .geojson and build search filter rather than set a static geojson_geometry variable. However, it did not work as it used to with the static variable. Win7 pro, Py 2.7.15 32bit. Below are the codes and error.

aoi = r"E:\PLANET_Py\multipart84.geojson"

with open(aoi) as f:
    data = json.load(f)

for feature in data['features']:
    poly_type = feature['geometry']['type']
    coords = feature['geometry']['coordinates']
    geojson_geometry=json.dumps((feature['geometry']))

geometry_filter = {
  "type": "GeometryFilter",
  "field_name": "geometry",
  "config": geojson_geometry
}

t = dt.now()

y = t - timedelta(days=1)
yy = t - timedelta(days=7)

# get images acquired within a date range from current date
date_range_filter = {
  "type": "DateRangeFilter",
  "field_name": "acquired",
  "config": {
    "gte": "%sT00:00:00.000Z" % (yy.strftime('%Y-%m-%d')),
    "lte": "%sT00:00:00.000Z" % (y.strftime('%Y-%m-%d'))
  }
}

# only get images which have <50% cloud coverage
cloud_cover_filter = {
  "type": "RangeFilter",
  "field_name": "cloud_cover",
  "config": {
    "lte": 0.5
  }
}

# combine our geo, date, cloud filters
combined_filter = {
  "type": "AndFilter",
  "config": [geometry_filter, date_range_filter, cloud_cover_filter]
}

search_result = requests.post(
    'https://api.planet.com/data/v1/quick-search',
    auth=HTTPBasicAuth('myAPIkey', ''),
    json=search_request)

print(json.dumps(search_result.json(), indent=1))

AND the error is:

>>> print(json.dumps(search_result.json(), indent=1))
{
 "field": {
  "filter.config.0.config": [
   {
    "message": "u'{\"type\": \"MultiPolygon\", \"coordinates\": [[[[176.46349778345987, -38.11473809258582], [176.46250239004382, -38.11601884487258], [176.4572092610443, -38.112539947133165], [176.45354029872024, -38.10990685811454], [176.45037426916238, -38.10994039692455], [176.4481870229376, -38.108336173676946], [176.44620341737033, -38.106249209991965], [176.4505735765721, -38.10296254706899], [176.46137920964145, -38.11140400374731], [176.46349778345987, -38.11473809258582]]], [[[176.45454303454335, -38.08318117683172], [176.45779400238163, -38.091368419149184], [176.45855077571625, -38.09462355135948], [176.45431160331216, -38.09427167466392], [176.4466842623444, -38.095628113517066], [176.4462098387771, -38.09355640167274], [176.4526823046912, -38.09169766808866], [176.45147933287407, -38.08696586891919], [176.4451995665532, -38.08810379555714], [176.44410615101478, -38.0856926990337], [176.45454303454335, -38.08318117683172]]]]}' is not of type 'object'"
   }
  ]
 },
 "general": []
}

I suspect that the unicode in python is problem

"u'{\"type\": \"MultiPolygon\", \"coordinates\":..........

But do not know how to avoid it. I thought using json.dumps() was enough.


texts: QGIS Multiple CASE WHEN/THEN Statements for Expression Based Labels

texts:
I am trying to label a single point with multiple labels, based on the attributes for that point.  For example, if one of the attributes is equal to 1, I want to label it with a specific text string; and if another attribute is equal to 0, I want to label it with a different text string, etc.  In addition, I want to use other attribute information to control the placement of the labels around the point, as well as their color.

Can multiple CASE WHEN/THEN expressions be stacked, to achieve the desired result?  I have tried to include multiple WHEN/THEN statements inside my CASE expression, looking like this:

CASE
WHEN attribute1 = 1 THEN 'Apple'
WHEN attribute2 = 0 THEN 'Grape'
WHEN attribute3 >= 0 THEN 'Pear'
END

This returns the desired label for attribute1, but ignores the labels for the other two attributes.  I believe this is because a CASE statement can only refer to a single variable/attribute at a time?  If that is true, is there another conditional expression type I can use to achieve what I want?  I'm trying to avoid having to add the point layer multiple times to my project, in order to label based on different attributes.


texts:
I have to add a point on a existing spplot in R. I need it only for a visual purpose. I have the coordinates of the point, but I don't know how to change the code in order to add the point.

This is the code of the spplot:

spplot(pioggiaTP[1])

And this is the map obtained by running the command:

The coordinates of the point are X=720731,8704 and Y=4861300,795 in the reference system WGS 84 / UTM zone 32N EPSG:32632, the same system of the plot.

Anyone can help me?


texts: Oracle 11g aliases that work in ArcMap 10.4

texts:
I am getting this error shown in the picture when I try to export a database view into a GDB table in ArcMap 10.4. I know what is causing the error, its the double quotes I use to create alias names for the columns that I am trying to export in the view. The reason I am doing this is because the company policy states that columns in datasets that we upload to ArcGIS Online must be plain text with no underscores, or other characters. So what im wondering is whether there is a way to create the aliases in Oracle 11g, PL-SQL Developer, that will work in ArcMap. Single quotes don't work they are an invalid identifier in Oracle 11g.

Here is an example of what is not working in ArcMap:

t.AS_BUILT_AGGREGATE_SIZE as "AS BUILT AGGREGATE SIZE",
t.as_built_hamburg_voids as "AS BUILT HAMBURG VOIDS",
t.as_built_vma as "AS BUILT VMA",
t.as_built_vfa as "AS BUILT VFA",

When I remove the alias, it exports fine. There are many columns and I would rather not do this by hand in ArcMap.


texts: Update query with sub query and join in Spatialite

texts:
I am looking to do an update query to one table (munparcel_zoning) with the sum of values from a column in another table (permability.permarea) where the type (permtype) is equal to a value (impervious). As this is a many-to-one, I am thinking I need a group by (parcelid). As there are two tables involved I am thinking the best way is to do an inner join, using parcelid = parcelid from each table. However, I am having issues with the query executing.

UPDATE munparcel_zoning
SET imperv =
(SELECT Sum(permability.permarea)
    FROM permability
WHERE permability.permtype = 'impervious'
GROUP BY permability.parcelid)
INNER JOIN munparcel_zoning
ON permability.parcelid = munparcel_zoning.parcelid;


texts: Feature Templates disappear in Arcmap after i remove Feature Class from Table of Contents

texts:
I am working on some Polygon Data in Arcmap and I am trying to generate about 20 or more feature templates for each Feature Class in the Geodatabase. All my feature classes are in one Feature Dataset and subsequently in a mother Geodatabase. I am able to create feature templates for each one and once i drag the Feature Class into my working area I see all of its associated feature templates in the 'Create Features' pane.

However once, I remove the particular feature class from Table of Contents and drag it back later on into the Working area from Catalog, all my generated feature templates vanish. Even when i restart the arcmap, I return to see all my created feature templates disappear.

Please help me out.


texts:
Possible Duplicate:
  GIS Apps for Android Tablets

Do you know any GIS applications for android OS?

I can recommend following ones:

Locus - it is equivalent of GPS tourist navigation

Shape file viewer - it is .shp .dbf reader

Arc Gis

It is also seems that Quantum GIS is starting to test for some devices.


texts: GIS Layer of USA State Court Jurisdiction Boundaries

texts:
Does anyone know of a national layer?

Edit: Sorry for the delay in clarifying.

I am looking for all courts (district, circuit, municipal, Superior, etc.) depending on the terminology used for the state for the USA.

Based on the answers provided and research I've done I don't believe this is available.

I very much appreciate all the feedback and leads.


texts: While populating table of memory layer one record is missing

texts:
I have script that calculates length of lines and adds values to attribute table. It works for vector layer, but for memory layer it skip always one value, as you can see on the pic. I calculated length by field calculator to show that it moves the value one row up.

Here  is the code:

from PyQt4.QtCore import *
layer = qgis.utils.iface.activeLayer()
provider = layer.dataProvider()
lengths = [feat.geometry().length()
        for feat in layer.getFeatures() ]

field = QgsField("length5", QVariant.Double)
provider.addAttributes([field])
layer.updateFields()

idx = layer.fieldNameIndex('length5')
for length in lengths:
    new_values = {idx : float(length)}
    provider.changeAttributeValues({lengths.index(length):new_values})

How should I edit code to get right value in right row? I suppose that the problem is in last line: the "{lengths.index(length):" doesn't work proper for memory layer. Am I right?

(Working in QGIS 2.12)


texts: FME Data Tranform and Concatenate Tables in DB without Key from Access to SQL

texts:
I have been able to do a conversion (in FME2013) from an Access DB to a SQL Express DB(in Desktop Workbench), adding some transformers to grab data from three of the tables, two of which have primary keys, and running a SQL query executer I can populate a new table;  querying out needed information from all three, (or inner join on the two with the key) this is all done one DBs at a time because the primary keys repeat.

I would now like to grab multiple Access DB’s all at once and do the same thing  only grouping them together based on the DBs they originate from adding a new primary key field containing the file name of the Access DB in that field.

When I use the Full Filename Attribute from the FilemnamePartExtractor it combines all of the selected filenames adding special characters (“,**) I have removed them using the StringReplacer.
I am at a point where I cannot seem to move forward successfully, due to my limited experience working with FME.

I am also wondering if this is the correct approach or is there a more straightforward method?

In the first and third writer I am using the FilemnamePartExtractor, StringReplacer, StringConcatenator and the BulkAttributeRenamer and dropping the values into a new SourceData… field.


texts: I'm building a municipal model with many pipes, manholes & other structures as vector layers - these are coming from several different sources of information and the field structure can be slightly different depending on the source.  I then need to add many new fields to each table for field data capture.  What I'm wondering, is if there's a way to build an empty table within a template, or a script that I can run (or something else) to automate this task?  I'll have many tables that I need to add these same fields to in this project, and I've got another dozen similar projects after this one... I would really like to find a way to add these same fields automatically rather than manually.


texts: Programmatically changing line lengths in QGIS?

texts:
Provided a set of lines with various lengths, I would like to change their lengths (hence, changing their geometry) to a specific value, using their centroid as base point.

I have checked the algorithms and whereas there is a similar one called Rotate, I haven't seen anything like what I am looking for.

Is there any way to achieve this?

I assume I could do it with python console, but I do not have a clue on what to type in order to achieve it - I am not familiar with python yet.


texts: ArcMap 10.3 returns DIST_2D in Skline Graphs as angular distance, not linear distance

texts: The Skyline processing frequently fails with error 999999 and others (050016, I believe).  When I can produce a Skyline feature and then produce a Skyline graph, the horizontal distance returned (titled DIST_2D in the table) is an angular distance, not a linear distance.  The elevation is therefore calculated as an arctangent of the linear vertical distance over the angular horizontal distance, a truly "interesting" result!


texts: AttributeError: 'NoneType' object has no attribute 'bounds'


texts:
I'm trying to read a kml from this kmz I unzipped
however, when I try to read a layer from it (as I've done it with previous kmzs) it says:

library(rgdal)
unzip("EstacionesClimatologicas.kmz")
ogr <- readOGR("doc.kml", "Puebla")

which gives the following output and error:

OGR data source with driver: KML
Source: "doc.kml", layer: "Puebla"
with 117 features
It has 2 fields
Error in readOGR("doc.kml", "Puebla") : Feature mismatch

when I check

  ogrInfo("doc.kml", "Puebla")

says it's wkbPoint with 3 dimensions,
How can I get this to run?


texts: Selecting maximum value in attribute table using Python code?

texts:
I saw a lot of posts about how you select maximum value in the particular field and I've been trying to use that query in a python code, but none of the posts really helped me solve my problem. There is a problem that I always get an error saying that I have an invalid expression.

I've been using:

arcpy.Select_analysis (input_layer, output_layer,"Shape_Area = (SELECT MAX( Shape_Area)" FROM layer_name)

I'm so confused as this query sometimes works and sometimes not, even in the attribute table or if I'm trying to use a tool Select, Table Select etc.
Can anyone help me? Or is it there some other way how to identify the maximum value in the field.


texts: Converting ArcMap Route Elements to KML?

texts:
I have used ArcMap to create bus routes with stops using the 'Find Route' tool.

I have saved my maps as MXD files and need to convert to KML.

I tried to convert to KML using the conversion tool but it didn't work.


texts:
I am using ArcGIS 10.2.2 and I need a way to automatically delete parts of a sketch. For example, I have multiple polygons with numbers inside of them. As shown below, the polygon has the number 28, and under the edit sketch properties window, the number 2 is listed under part 1 and the 8 is listed under part 2. I feel like there should be a way to select multiple polygons, and have a script that could delete any vertices that are not included in part 0, that way I would not have to double click each polygon to delete superfluous vertices.

Can this be solved via a tool in ArcGIS?

I want to remove the numbers, not generalize the border.


texts: How can I access WMS service from Geoserver using ArcGIS Silverlight API?

texts:
Having GeoServer installed locally in my machine, I'm currently doing a R&D by trying to access WMS service from GeoServer using ArcGIS API for Silverlight.

I'm trying to replicate the sample in ESRI site which deals with accessing a WMS service using ArcGIS API for Silverlight. I'm unable to view the map in the browser.

Given below is the XAML code used in my application:

<UserControl x:Class="SilverlightApp.Beginners.MainPage"
  xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
  xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
  xmlns:esri="http://schemas.esri.com/arcgis/client/2009">
    <Grid x:Name="LayoutRoot">
        <esri:Map WrapAround="True">
            <esri:WmsLayer
              ID="OtherWMSLayer"
              Layers="OSOpenData"
              Url="http://localhost:9595/geoserver/wms"
              Opacity="0.7"
              SupportedSpatialReferenceIDs="4326"
              SkipGetCapabilities="True"
              Version="1.1" />
        </esri:Map>
    </Grid>
</UserControl>

Please can anyone shed light on this?


texts: constructing customer-derived trade areas

texts:
I am interested in creating some customer-derived trade areas that represent "most" of the customers for a given store (say 80%). I have both the customers and the stores geocoded.

The ESRI Business Analyst description of the general idea (with pictures) is here. Unfortunately I work at a MapInfo shop, so I am trying to hack it together myself. Does anyone know if something like this exists already so that I don't have to reinvent the wheel?

Here's a summary of what I've tried so far:

I've start at the store, ranking the customers from closest to farthest, and then adding each customer in until I get to ~80% of sales for each store. Then I put the convex hull around those folks and call that the trade area. This works reasonably well, but it misses major customers who happen to be distant.

I've also tried ranking the customers in terms of sales, from largest to smallest. Then I keep the ones who make up the 80%, and then snap a convex hull around them. This gets the big guys, but I loose the 80% feature since that area now includes lots of the small shrimp.

I've also tried a hybrid approach where calculate z=sales/(straight line distance) for each customer, rank by z, and put CH around the folks that add up to 50%. That CH will also include other customers that were not selected, and so I then trim the farthest customers in the hull until I get back down to 80% of sales. This works the best of the 3 methods.

Any advice would be greatly appreciated.


texts: Python window does not open in ArcGIS 10.4

texts:
I have a Python script that I need to run in the Python window. When I click on the Python window icon, nothing happens. I expect a window to open up which allows Python commands to be entered. I think this may be due to removing additional monitors. Maybe the Python window is opening on a phantom monitor?

Is there a way for me to get the window back on the main monitor?


texts: Speeding up length calculation in QGIS 2.0 Python plugin

texts:
In a qgis 2.0 plugin for windows, I calculate length and update a field with the values. It works, but takes a lot of time. Execution for this part only takes about ten minutes. Can anyone suggest anything to speed things up?

def calculatelength(self):
    idx = self.resLayer.fieldNameIndex('Longueur')
    caps = self.resLayer.dataProvider().capabilities()
    features = self.resLayer.getFeatures()
    for feature in features:
        if caps & QgsVectorDataProvider.ChangeAttributeValues:
            fid = feature.id()
            flen = feature.geometry().length()
            attrs = {idx : flen}
            self.resLayer.dataProvider().changeAttributeValues({fid : attrs})


texts: Merging/Dissolving polygons by Common Attribute Field in QGIS?

texts:
I have joined corresponding Local Enterprise Partnerships table data into a polygon shapefile containing UK Local Authorities.

How do I Dissolve or Merge the the Local Authorities polygons by LEP_NAME as shown in the screenshot?

Ideally, I would like to retain the Local Authority attributes with the resultant merged polygons. I have zero scripting knowledge so am looking for a tool based solution. Any help or pointers much appreciated.

Additional Note: When attempting to use the dissolve tool the LEP_NAME field is not available as shown above. The LEP attribute data was created in Excel, saved as CSV format and joined to the Local Authority polygon shapefile. There are duplicate IDENTIFI0 fields as shown in the other screenshot - could this be the root of the issue?


texts: Create polygons from many small. Expand to include space between

texts:
I'm just a programmer without GIS-skills as such. Hence this question:

Intro:

I have a lot of polygons; actually each private-owned real estate in my country.
I also have "counties" as polygons.
In each county there are a varying number of (administrative) subdivisions (2 to 55) - and each real estate belongs in one of these subdivisions.
But real estates are divided by roads and there are blank areas without any relation to the subdivisions.

My goal is

To have just one polygon for each subdivision.
Have the blank spaces between real-estate-polygons divided evenly between subdivisions.
The county marks the outer border.

Simple example:

All the red boxes are merged. All the green are merged. The spaces between are shared. The black box is the outer limit.

Is there a readily available tool that can do this easily? Like QGis, ArcGIS, a db-engine or library/API?
I've tried perusing the different options, I've tried testing various approaches, but to no avail.

If there aren't any tools, I guess I need to program it myself. Any thoughts?

If there is a term describing the process above, please let me know: I'm such a GIS-noob, so I don't even know what to Google :-)


texts:
I am looking for training app installed inside QGIS, which can be used for training purposes.

I mean, a plugin like an e-learning software than can drives you directly inside the QGIS software and can check your result automatically.

I looked under python plugin repository with tutorial or training keys but with no good result.

I think it should be better integrate with a C++ plugin or a specific QGIS standalone App than a python plugin, but I could be wrong.

I am also interest about project which use a kind of tutorial module, like step by step, with highlighting zone, mouse move, button action, setting parameters directly under QGIS.

If I am still not clear enough, this kind of app are usual on web, for e.g., when you launch openstreetmap editor on their website, there is a kind of "beginner tutorial". I would know if someone develop something similar under QGIS desktop, never mind if it's on a totally different way.
it´s a kind of state of art about tutorial existing tools under QGIS.


texts: qgis - slope from DEM - errors due to projection and oversimplified 'z factor'

texts:
I'm trying to generate a slope raster from a DEM raster.  It seems there is a pretty significant error due to the 'z-factor' option.

I'm fairly new to this, but some other posts and responses indicate a general workflow:

1) obtain a DEM.  The ones I got (NED) are unprojected: horizontal coords are lat-lon, and vertical coords are meters.

Two possible branches here, A or B:

A2) project the DEM to UTM so that horizontal and vertical coords are in the same units
A3) then do slope analysis on the UTM-projected DEM, using z-factor = 1.0

OR

B2) do slope analysis on the unprojected DEM, using z-factor = 111120 (which has been suggested in some posts as a number of meters per degree - see error discussion below)
B3) then project the raster to UTM

Problem with option A: the projection process introduces gridded artifacts, and these get amplified in the slope raster:

Original DEM (NAD83, unprojected; zoomed in to one peak; min/max set to show high contrast here):

DEM projected to NAD83 UTM zone 10, note gridded artifacts:

slope analysis output from NAD83-UTMz10-projected DEM, note the same artifacts:

Problem with option B: 111120 seems vastly oversimplified and produces significant errors, see error table below - though, gridded artifacts are greatly reduced since they don't exist in the DEM that you run slope analysis on.  First, number of meters per degree of latitude is not the same at all latitudes, due to the earth not being a perfect sphere.  Second, and more obvious, number of meters per degree of longitude is a function of latitude, varying from ~111km at the equator to 0 at the poles.

For our area (Nevada County, California) here are some numbers from an online calculator, referenced in other places as a trusted source (but feel free to verify these):

Nevada County S extent ~ 39.0deg (WGS84)
Nevada County N extent ~ 39.5deg (WGS84)
Nevada City ~ 39.25deg (WGS84)

39.0 deg:
1 deg lat = 111015.45 meters
1 deg lon = 86626.37 meters

39.25 deg:
1 deg lat = 111020.23 meters
1 deg lon = 86320.71 meters

39.5 deg:
1 deg lat = 111025.01 meters  (about 9.5 meters longer than at S extent of Nevada County)
1 deg lon = 86013.39 meters (about 613 meters shorter than at S extent of Nevada County)

(at equator, for reference:)
0.0 deg:
1 deg lat = 110574.27 meters
1 deg lon = 111319.46 meters

So, making some straightforward assumptions about how qgis uses the 'Z Factor' argument, here are some resulting errors (note this is in Nevada County but not the same location as the above images), along with Range-Bearing Line from Terrain Navigator Pro which is treated as gospel in this case; also note that this is just for one sample set of from/to; error would vary by compass direction of the from-to line also:

(z-factor is entered by hand for each row; all the other columns are calculated from it)

So, any ideas?  Am I calculating something wrong, or is there a better workflow, or is there just no really good way to get around the error?


texts: How to return array from label offset expression?

texts:
I am trying to set the x,y offset of layer labels based on a data defined override. What I'm trying is the following:

CASE
WHEN name is 'Brooklyn Bridge' [3.5, 0]
ELSE [-3.5, 0]
END

However this does not seem to be the correct syntax to set x and y values simultaneously. How would you do this?


texts: Add Data Widget - prevent added layers from showing legend?

texts:
Is there a way to prevent a "curated" web layer from being automatically included in the legend after it has been added to the web app with the Add Data Widget?

Background:

We have a lot of layers on our Web GIS and decided to de-clutter portal maps/apps Layer Lists by removing layers that aren't frequently used but still allow users to add them into the web apps by using the Add Data Widget (after marking select layer services as "curated")

The Add Data Widget works well and users can add any of the "curated" layers, however these added layers get automatically added to the legend although we've configured the layers to hide legends and saved the layers (Save Layer) in portal. I was expecting the same behaviour as if I were to add a layer that has been configured and saved in portal to a new map - a previously saved layer inherits the hide in legend setting, popups and other behaviour. It seems when such layer is added to a RUN-TIME app via the Add Data Widget it ignores the the saved layer behaviour: shows the legend by default and ignores any custom popup behaviour.


texts: Multiband virtual raster with gdalbuildvrt: Control how bands are sorted

texts:
I am making a multiband virtual raster using gdalbuildvrt. Everything works fine, but a little detail remains annoying: I cannot control how the bands are sorted. I specified the input raster using the -input_file_list option, but the order is not kept by gdalbuildvrt. Instead it sort the band with lexicographic sorting, which is not what I want.

For instance, that would be my input file list:

/home/example/rasters/1634.tif
/home/example/rasters/5899.tif
/home/example/rasters/15680.tif

I'd  like to have a vrt corresponding to:

band 1 -> 1634.tif
band 2 -> 5899.tif
band 3 -> 15680.tif

But instead I get:

band 1 -> 15680.tif
band 2 -> 1634.tif
band 3 -> 5899.tif

I could manually rearange things, but I'd like to avoid it if there is any cleaner solution...

So my question is: is there a way to control how gdalbuildvrt sort the bands when creating a multiband virtual raster?


texts:
I use FME a lot for manipulation of spatial data, and would like to leverage it's Python library, fmeobjects in PyQgis.

If I run the following in my standard Python IDE it works fine:

import sys
sys.path.append("C:\\Program Files (x86)\\fme\\fmeobjects\\python27")
import fmeobjects

But the exact same code, when run in PyQgis, throws
    ImportError: DLL load failed: The specified module could not be found.

In actual fact, fmeobjects is a .pyd file. I don't know if that's the problem, as my understanding is that pyd and dll are analogous to each other.

How can I get fmeobjects to import into PyQgis?


texts: Creating planar 3D buffer from 3D polyline feature class in ArcGIS Desktop

texts:
I want to create a planar (=without volume) buffer polygon or a TIN from a polyline feature class that represents a road centreline ,containing Z-value. I want to expand my polyline horizontally (=parallel to ground) in 3D environment to get a 3D polygon which its volume is zero. I tried to use buffer analysis tool and enabled z value, but it results in a polygon without elevation. What I want is as same as Creating a 3D buffer polygons (plane) from a 3D polyline that follow line catenary curvatures except my goal is to calculate the cut & fill for a road.

In the attached image, the black polyline represents road centreline. The black polygon is the result of buffering and the red polygon, is a polygon which represents terrain. My ultimate goal is to calculate the volume between two polygons or surfaces to get cut & fill.

How can I create my desired polygon or surface in ArcGIS Desktop?

I have a 3D Analyst license and am using ArcMap/ArcScene/etc.


texts: Updating file output for model based on parse path of raster iterations

texts: I have created a model in ModelBuilder that uses the iterate tool (recursive) to search through a file directory and locate rasters for input into the model. Where my problem arises is in the output at the end of the model. I am trying to have the model output the file in the sub folder from wich the iteration pulled the raster. Below I have included images of my ModelBuilder.


texts: Performing Random Forest using EnMap Box?

texts:
I am new to the software EnMap Box that incorporates Random Forest (RF) as one of its capability.  To learn the Enmap, I followed the RF tutorial provided online.

My goal is to use RF as a tool to classify Landsat images using more than 100 explanatory variables.

Link to RF: https://bitbucket.org/enmapbox/enmap-box/wiki/imageRF%20-%20Manual%20for%20Application#!variable-importance-of-rfcrfr-models

In the EnMap tutorial (link provided above), under the "Parameterization of RFC/RFR Models", it only says "Select the Image to be classified...." But what image in particular? Am I going to use the image that stacked all 100 variables that I want in the analysis?  In other words, how would I input the various explanatory variables into EnMap, so that I could then apply the tool "Variable Importance" to select important layers for classification?

Is there anyone familiar with EnMap and could guide me through this?


texts: Uploading fields from Excel into ArcMap Attribute table?

texts:
I have digitized a polygon using ArcMap 10.1.

I want to upload fields from an excel table into the attribute table of the polygon layer.

How do I do this ?


texts:
I have experience in PostgreSQL but I am fresh new in PostGIS and I can see how quickly would I like it. I have extensively worked on time series and now I am discovering the spatial side of my data. I have installed PostGIS on my server and then I am wondering how to anwser tons of questions. Thus I am in the prospective part of my work and I have everything to learn.

Anyway, I would like to have some clues from GIS expert in order to take the good start. Basically, I am monitoring concentrations at constant places, but we also have some campaigns where we monitor concentrations when we are in motion. That is, we have, at least, a five dimension problem, spatial and time is tracked down by GPS and another dimensions are assessed by dedicated monitors.

I know that PostGIS have a multipoint object, but all example that I have seen on the web are limited to 4D. Is it possible to use this structure to store multiple concentrations in addition of space and time. If so, is there a simple way to get rid of extra dimensions such as concentrations and time, when we are just working on space? If I want to post process my data with QGIS, will I be able to use such structure?


texts: How to add command developed with C# to ArcToolbox?

texts:
I have a command that I developed with c#.

I want to add this command to the ArcToolbox.

Is there any way to add such a command to ArcToolbox in ArcGIS Desktop programmatically ?


texts: postgisQueryBuilder in QGIS can't connect to PostgreSQL DB - QPSQL driver not loaded

texts:
PostgreSQL is installed properly and I can work with my DB very well. My databases show up as connected inside the QGIS browser window. I think the appropriate packages are installed:

# rpm -qa post* libqt*
postgresql9.6-contrib-9.6.15-1.mga6
postgresql9.6-9.6.15-1.mga6
postgis-2.2.2-3.mga6
postgresql-jdbc-9.4.1210-1.mga6
libqt4-database-plugin-sqlite-4.8.7-15.1.mga6
postgresql9.6-plpgsql-9.6.15-1.mga6
libqtsql4-4.8.7-15.1.mga6
postgresql9.6-server-9.6.15-1.mga6
libqt4-database-plugin-pgsql-4.8.7-15.1.mga6

Can someone interpret this stacktrace for me? This is what I get when I try to connect to my DB inside postgisQueryBuilder:

Warning: QSqlDatabase: QPSQL driver not loaded
Stacktrace (piped through c++filt):
src/core/qgsmessagelog.cpp: 45: (logMessage) [9835ms] 2019-11-14T17:25:34 Qt[1] QSqlDatabase: QPSQL driver not loaded
qgis[0x40b7c4]
qgis[0x40bafa]
/lib64/libQtCore.so.4(qt_message_output(QtMsgType, char const*)+0x16)[0x7f46c3634356]
/lib64/libQtCore.so.4(+0x6f529)[0x7f46c3634529]
/lib64/libQtCore.so.4(qWarning(char const*, ...)+0x94)[0x7f46c3634774]
/lib64/libQtSql.so.4(+0x14a27)[0x7f46c1ab5a27]
/lib64/libQtSql.so.4(QSqlDatabase::addDatabase(QString const&, QString const&)+0x11)[0x7f46c1ab5d41]
/usr/lib64/python2.7/site-packages/PyQt4/QtSql.so(+0x47c8d)[0x7f4640c65c8d]
/lib64/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8a84)[0x7f4643d039b4]
/lib64/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x667c)[0x7f4643d015ac]
/lib64/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x85c)[0x7f4643d048dc]
/lib64/libpython2.7.so.1.0(+0x70a6c)[0x7f4643c8ea6c]
/lib64/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7f4643c6a8b3]
/lib64/libpython2.7.so.1.0(+0x5af47)[0x7f4643c78f47]
/lib64/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7f4643c6a8b3]
/lib64/libpython2.7.so.1.0(PyEval_CallObjectWithKeywords+0x47)[0x7f4643cfa927]
/usr/lib64/python2.7/site-packages/sip.so(+0x13738)[0x7f464380d738]
/usr/lib64/python2.7/site-packages/PyQt4/QtCore.so(PyQtProxy::invokeSlot(qpycore_slot const&, void**, int)+0x153)[0x7f464351de43]
/usr/lib64/python2.7/site-packages/PyQt4/QtCore.so(PyQtProxy::unislot(void**)+0x105)[0x7f464351dc1f]
/usr/lib64/python2.7/site-packages/PyQt4/QtCore.so(PyQtProxy::qt_metacall(QMetaObject::Call, int, void**)+0x8a)[0x7f464351db02]
Warning: QSqlDatabase: available drivers: QSQLITE
Stacktrace (piped through c++filt):
src/core/qgsmessagelog.cpp: 45: (logMessage) [6ms] 2019-11-14T17:25:34 Qt[1] QSqlDatabase: available drivers: QSQLITE
Driver not loaded Driver not loaded
qgis[0x40b7c4]
qgis[0x40bafa]
/lib64/libQtCore.so.4(qt_message_output(QtMsgType, char const*)+0x16)[0x7f46c3634356]
/lib64/libQtCore.so.4(+0x6f529)[0x7f46c3634529]
/lib64/libQtCore.so.4(qWarning(char const*, ...)+0x94)[0x7f46c3634774]
/lib64/libQtSql.so.4(+0x14aad)[0x7f46c1ab5aad]
/lib64/libQtSql.so.4(QSqlDatabase::addDatabase(QString const&, QString const&)+0x11)[0x7f46c1ab5d41]
/usr/lib64/python2.7/site-packages/PyQt4/QtSql.so(+0x47c8d)[0x7f4640c65c8d]
/lib64/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8a84)[0x7f4643d039b4]
/lib64/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x667c)[0x7f4643d015ac]
/lib64/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x85c)[0x7f4643d048dc]
/lib64/libpython2.7.so.1.0(+0x70a6c)[0x7f4643c8ea6c]
/lib64/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7f4643c6a8b3]
/lib64/libpython2.7.so.1.0(+0x5af47)[0x7f4643c78f47]
/lib64/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7f4643c6a8b3]
/lib64/libpython2.7.so.1.0(PyEval_CallObjectWithKeywords+0x47)[0x7f4643cfa927]
/usr/lib64/python2.7/site-packages/sip.so(+0x13738)[0x7f464380d738]
/usr/lib64/python2.7/site-packages/PyQt4/QtCore.so(PyQtProxy::invokeSlot(qpycore_slot const&, void**, int)+0x153)[0x7f464351de43]
/usr/lib64/python2.7/site-packages/PyQt4/QtCore.so(PyQtProxy::unislot(void**)+0x105)[0x7f464351dc1f]
/usr/lib64/python2.7/site-packages/PyQt4/QtCore.so(PyQtProxy::qt_metacall(QMetaObject::Call, int, void**)+0x8a)[0x7f464351db02]

Do I somehow have to load that QPSQL driver before? Or do you see hints about missing packages inside this? I'm relatively new to this.


texts: How to convert GeoTIFF to grayscale and add gaussian blur?

texts: I have a giant GeoTIFF I'd like to convert to grayscale and add a gaussian blur to. As far as I know, GDAL can strip two of the three color bands and call it grayscale, which might work, but is there a better way? Imagemagick has both grayscale and gaussian blur functionality, but can not process georeferenced TIFFs (as far as I know). Any suggestions how I should go about doing this?


texts: forEachFeatureAtPixel not limiting features?

texts:
I'm trying to limit the features returned when I click on the map to just those in the endPoint layer:

    endPoint = new ol.layer.Vector({
            source: new ol.source.Vector({
                features: endPoints //add an array of features
            }),
        });
    map.addLayer(endPoint)

Unfortunately, before this layer is even initialised, whilst clicking on features in other layers they are getting logged to the console. I feel like I'm using the correct syntax to limit the features to just those in the endPoint layer (forEachFeatureAtPixel(pixel, callback, opt_this, opt_layerFilter, opt_this2)) but I must be doing something wrong - if I try to return layer rather than feature it gives me a null value?

Is the issue that I need to name the layer somehow?

    map.on('click', function(evt) {
      var pixel = map.getEventPixel(evt.originalEvent);
      var mapFeature = map.forEachFeatureAtPixel(pixel, function(feature, layer) {
          return feature;
            }, null, function(layer) {
                     return layer === endPoint;
            });
      console.log(mapFeature)
      })


texts: Exporting DEM-derived rasters to CSV from QGIS?

texts:
I have been able to import into QGIS a set of raster topographic maps and the associated DEM raster maps. It registered the two for me. From there I have calculated the associated slope raster map. Now I want to export the data to a CSV file (or similar format) for processing in a program that I am writing. Ideally the data would be structured in a way that would permit me to query for the slope and elevation data for any given point on the map.

If this is possible, then how do I do it?

Specifically, I have been able to export the data (slope & DEM) to two different CSV files. However I would prefer the data in one file if possible.


texts: QGIS Download getting stuck, NSIS error?

texts:
I need to edit a DBF file, and I tried to install QGIS at http://www.qgis.org/en/site/forusers/download.html.

However, none of the standalone installers work. I tried installing 32/64bit, Latest/Long-term release, it didn't matter.  It would get stuck downloading the .exe file itself for about a half hour before it would seemingly finish. But then when I tried to run the file, get:

Then I found this and I tried to install it using the OSGeo4W Network Installer, but even that got stuck at 6%. Is there anything I can do to solve this issue?? I tried this on two different computers and got the same result...


texts: QGIS : See ID in polygons in ascending order

texts:
I have input polygons in shp, I assign id with $rownum but to display a label with the order isn't good. It seems to bind to the entry order.

How to automatically change the id for them in order to display?


texts: OpenLayers Popup relative to a feature's position - prevent overlaps - only one popup in one time

texts:
There are many questions related about this topic. However there is no complete example for this widely used pop-up type:

Basically i need a popup which will appear the top-right of a Feature (Point) when user clicked on that feature. When i am able to do this, unfortunately pop-up overlaps part of the feature. So probably i should set an offset to the feature.

Additionally when map is clicked, the popup should be removed and also when another feature is clicked, the old popup should be removed as well as new popup has been created.

Which event listeners should i use FeatureSelected or FeatureClicked?

Which popup type should i use Anchored or something else?

Could this be a projection issues problem?

I just started to use openlayers and javascript.


texts: Adding census tracts to addresses using ArcGIS for Desktop?

texts:
I am as rookie as they come in this space. I have a large list of addresses and longitude latitude points for different locations within a state. What is the easiest and most efficient way to determine the respective census tract number for each address? Over 50,000 addresses and I just want to pull in the tract into excel for each

Right now it is just one single data set. I can manually look up the address in some third party systems, but that process is extremely inefficient. I also don't know how to use arc gis but I've tried valiantly to import the data points into the map, then layer the census tract map over that. Unfortunately, I don't know if there is a way to merge those two sets of data then export from ArcGIS back into the spreadsheet?


texts:
I'm trying to compile OSRM in a Xubuntu 12.04 virtual machine following these instructions. I had some problems installing boost (a requisite) and I compiled it manually.

The error is this:

/home/user/src/Project-OSRM/datastore.cpp: In function 'int main(int, const char**)':
/home/user/src/Project-OSRM/datastore.cpp:163:86: error: no matching function for call to 'min(long unsigned int, std::basic_string<char>::size_type)'
/home/user/src/Project-OSRM/datastore.cpp:163:86: note: candidates are:
/usr/include/c++/4.6/bits/stl_algobase.h:187:5: note: template<class _Tp> const _Tp& std::min(const _Tp&, const _Tp&)
/usr/include/c++/4.6/bits/stl_algobase.h:233:5: note: template<class _Tp, class _Compare> const _Tp& std::min(const _Tp&, const _Tp&, _Compare)
make[2]: *** [CMakeFiles/osrm-datastore.dir/datastore.cpp.o] Error 1
make[1]: *** [CMakeFiles/osrm-datastore.dir/all] Error 2

I've tried to define these variables in compile time

CFLAGS=/usr/local/include
LDFLAGS=/usr/local/lib

The directories where boost headers and libs can be found. Just in case. Still getting the same error.

Any suggestions? I've used boost 1.49 (I know there's a problem with it, but I'm still not there). And maybe the problem isn't boost, but something else.

Many thanks in advance


texts: Calculating nearest distance between points and Spatial Polygons in R?

texts:
So right now I have pretty large datasets, about 200,000 (lat,long) locations and a Large SpatialPolygonsDataFrame with a polygons list of length 428386 spread in a single city of around 450 square kms.

I want to calculate the nearest Spatial Polygon for every point in the dataset. I have tried the following procedures:

Using the gDistance function from rgeos package. First I converted both the points and Polygons to same CRS and to planar projection and applied the function. Currently, around 75 points take 1.25 minutes which comes down to around 55 hours for my dataset.
Using dist2Line function from geosphere package. This allows me to use my (lat, long) points as it is, i.e. without modifying the data but it takes a huge amount of time. A single point took around 15 minutes so it is clearly not feasible.

Can anyone suggest me a better way of calculating these distances, if possible in R?

However I'm also comfortable in using python if it has a better approach.


texts: Changing class breaks in ArcGIS without resetting whole symbology?

texts:
For the purposes of this map, I have a raster with very specific class breaks with custom colors for each class.

Sometimes for various reasons the class breaks have to be adjusted slightly based on new input.

However changing the classification breaks completely resets everything.

The colors and labels will all be reset and have to be redone manually.


texts: Zoom Canvas to Selected Features with Python call from QGIS Layer Action

texts:
I have been trying to get one layer's clicked-on feature to have an effect on another layer's viewed features with some success, but I am unable to get the map canvas to be filled by the the selected layer's features when using python code called from within the layer actions call.

For reference, I have two layers - one not referred to in the code is a layer of zones. The other (referred to below) is the same zones intersected with same to identify all zones with their neighbouring zones. My intention is that on clicking on a zone in the former layer, I get the neighbouring zones highlighted, and zoomed to as a seamless operation.

When I run this code from the python console, it behaves exactly as I was hoping for:

registry = QgsMapLayerRegistry.instance();
name= 'intersecting_zones';
layer = registry.mapLayersByName( name )[0];
layer.removeSelection();
expr = QgsExpression("\"msoa_code_\"='E02002869'");
it = layer.getFeatures(QgsFeatureRequest(expr));
ids = [i.id() for i in it];
layer.setSelectedFeatures(ids);
box = layer.boundingBoxOfSelected();
iface.mapCanvas().setExtent(box);
iface.mapCanvas().refresh();

This will not work in the exact form within the python call of the original layer's action properties and I have understood that instead of using iface, I needed to switch to using a reference to QgsMapCanvas, hence I use code of the following form - first easy to read, second same code in the single line for use within Layer Actions.

from qgis.coreimport QgsApplication;
from qgis.gui import QgsMapCanvas;
registry = QgsMapLayerRegistry.instance();
canvas = QgsMapCanvas();
name= 'intersecting_zones';
layer = registry.mapLayersByName( name )[0];
layer.removeSelection();
expr = QgsExpression("\"msoa_code_\"='[% "msoa_code_" %]'");
it = layer.getFeatures(QgsFeatureRequest(expr));
ids = [i.id() for i in it];
layer.setSelectedFeatures(ids);
canvas.setExtent(layer.extent());
canvas.refresh();

from qgis.core import QgsApplication; from qgis.gui import QgsMapCanvas; registry = QgsMapLayerRegistry.instance(); canvas = QgsMapCanvas(); name= 'intersecting_zones'; layer = registry.mapLayersByName( name )[0]; layer.removeSelection(); expr = QgsExpression("\"msoa_code_\"='[% "msoa_code_" %]'"); it = layer.getFeatures(QgsFeatureRequest(expr)); ids = [i.id() for i in it]; layer.setSelectedFeatures(ids); canvas.setExtent(layer.extent()); canvas.refresh();

This code (from within Layer Actions) works to a point, I get no error message and the desired target features are selected. However, I cannot get the zoom to selected features functionality to work at all. Am I using the pyqgis code incorrectly, or is there some other approach I should be using to get this last bit working?


texts: How to properly initialize ArcMap data frame created with IMaps.Create method?

texts:
In a custom ArcMap command implemented in .NET C# I create a new data frame this way:

// mxDocument is IMxDocument
IMap newMap = mxDocument.Maps.Create(); // create new data frame
newMap.Name = "I can name it as I want: " + Guid.NewGuid();
mxDocument.Maps.Add(newMap);
mxDocument.ActiveView = (IActiveView)newMap; // activate new data frame
mxDocument.UpdateContents();

However, compared to the data frame created with standard ArcMap means, there is something wrong with my data frame:

ArcMap doesn't show property page for the data frame (right click -
Properties)
It's not possible to remove my data frame from ArcMap
(right click - Remove)

Do I have to add some initialization code for IMap created with IMaps.Create (IMaps2.Create)?

The help for the method is quite laconic: "Creates a new map"...

I use ArcGIS version 10.2.2.


texts: How can I find the difference between two times (hh:mm) in QGIS?

texts: I have a shapefile of crime reports that includes the time a call was received and the time officers arrived on scene.  I want to create response times from the difference in these two values.  I know QGIS has support for the difference between dates, but I need the difference between two times in hh:mm format.  Is there a way to do this in QGIS specifically? I am running QGIS 2.0 standalone version.


texts: How to stack bands in Google Earth Engine?

texts:
I have created an Image collection in GEE and with the help of a function I have computed NDVI index and mapped it to create another collection with NDVI as a band.

Now I want the create a stacked image with NDVI bands of the entire image collection into one image. So it should be like NDVI_1, NDVI_2 and so on...

How can I do this? I am pasting the code that shows the NDVI collection that I have so far

// Collection of Images
var collection = ee.ImageCollection([feb1,feb2,Mar2,April1, April2, May1, May2, Jun1,Jun2,
July2, Aug2, Sep1, Sep2,Oct1, Oct2, Nov1, Nov2, Dec1, Dec2 ]);

//Using the following function,NDVI of the entire collection is computed
var indicesS2 = function(scene)
{ var ndvi = scene.normalizedDifference(['B8', 'B4']).rename('NDVI');
  var image = ee.Image()
                .set('system:time_start', ee.Date(scene.get('system:time_start')));
         return image.addBands([ndvi]).clip(Sheikhupura);
};
var NDVIcollection = collection.map(indicesS2);
print (NDVIcollection, 'NDVI');


texts: Copying representation symbology from one ArcSDE version to another?

texts:
We are having problems with our representations when copying from an ArcGIS 9.3 database to a 10.0. Representation symbologies with dash patterns (dots, dotted lines etc.) are reset to default settings in the version 10.0 database.

The problem is the same when you copy from a 10.x to 10.x database.

Have anyone experienced the same, and what did you do to make it go away?


texts: How to "page" through a set of features and display an infoWindow for each feature?

texts:
In the ArcGIS Server JS API the Popup object may be defined against an array of features.

In that case the Popup contains Next/Previous buttons to allow the user to "page" through each feature in turn, as in this example where there are multiple coincident features at the chosen location:

Is this possible using OpenLayers? That is, given an array of features can I easily navigate the infoWindow forward/backward to display the relevant information for each feature?

Thanks


texts:
I am working with arcmap 10.2.2.
I would like to label features in maps which I create with data driven pages. I have one ddp layer and another layer with the labels. The layer with the label features have map nummers (a field in the attribute table) according to the map numbers in the ddp. Now I only want those features to be labelled if the map number (on the dpp map) corresponds with the map nummer in the attribute table. Some maps have for example six line features, but only three should be labeled because in the attribute table only those three features have the corresponding map number. Now It labels all line features but the maps overlap so I do not want that.
Hope I made myself clear and someone knows How to do this. I tried already with the maplex label engine..

P.s. if possible NOT to difficult python scripting


texts: Can't get ESRI.ArcGIS.ADF.Connection.Local to work for 10.2.1

texts: Not sure what I'm missing. I'm trying to use the ESRI.ArcGIS.ADF.ComReleaser and I can't get C# to recognize the ADF namespace. I can add the reference by going to Reference->Add ArcGIS Reference, but when I build the project I get the error: "The type or namespace name 'ADF' does not exist in the namespace 'ESRI.ArcGIS' (are you missing an assembly reference?)". I saw another post that said you have to add ...ADF.Connection.Local as a .NET reference, but when I try that approach it's not visible in the list, only ...ADF.Local.


texts: Mosaic rasterstacks based on minimum of certain layer?

texts:
I am trying to mosaic 42 remote sensing rasterstacks (with 250 bands) based on the criterion that in overlapping areas, the pixel should be taken that has the most nadir viewing angle

Beside my rasterstacks I also have 42 rasters (so one for each stack) with the corresponding viewing angle for each pixel.

Any idea how to solve this?

I tried to include the viewing angle raster in the stack, and use something similar to

mosaic(a,b,fun=function(x)(min(x[[251]]))

but that didn't work...

Any advice?


texts: Printing max and min values of Landsat NDVI time series to GEE console

texts:
I am trying to extract the min and max values of an NDVI time series and print those values to the console, so I could then subtract those values to determine the seasonal stability of an area.  Most of the questions and answers in this forum explain how to reduce a collection to a Max NDVI composite image.  I just want the values.  I know I can download the CSV and get the values from there, but I am wishing to do this programmatically in Google Earth Engine.
Note: in the code below, the bounds are filtered by a polygon which I couldn't figure out how to attach.  Creating any sample polygon should work to run the code.
//This approach starts with a collection of L8 TOA images and filters them by cloud cover
  // Note: - CC is calculated only over the ROI (ie. not whole image CC score).
  //       - the simpleCloudScore algorithm only works with TOA images.

// --------------Set Location and Filter by Date and Cloud-------------------//

var MAX_LIST_SIZE = 100;

var L8_TOA = ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA")

    .filterDate('2013-12-25', '2016-12-25')
    .filter(ee.Filter.eq('WRS_PATH', 41))
    .filter(ee.Filter.eq('WRS_ROW', 36));
var count = L8_TOA.size();

var c = L8_TOA.filterBounds(geometry);

var withCloudiness = c.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: geometry,
    scale: 30,  });
  return image.set(cloudiness);});

// Here is where the computed cloud score threshold is set!
var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 5));

// print('L8 TOA Filtered Collection, by date and cloud',filteredCollection);

//---------------------------Add NDVI Band--------------------------------//

// Create an NDVI band
var addNDVI = function(filteredCollection) {
var nir = filteredCollection.select('B5');
var red = filteredCollection.select('B4');
var ndvi = filteredCollection.normalizedDifference(['B5', 'B4']).rename('NDVI');
return filteredCollection.addBands(ndvi);};

// Add the NDVI band to each image in the collection

var withNDVI = filteredCollection.map(addNDVI);

// print('L8 Filtered Collection with NDVI band',withNDVI)

// Filter collection to only include NDVI Band

var NDVI_bands = withNDVI.select('NDVI');

print('L8 Collection with cloud filtering added', NDVI_bands)

//--------------------------Create Map and Chart-----------------------------//

// Create an image time series chart, average it over the ROI

var chart = ui.Chart.image.series({
  imageCollection: NDVI_bands,
  region: geometry,
  reducer: ee.Reducer.mean(),
  scale: 200});

// Add the chart to the map.

chart.style().set({
  position: 'bottom-right',
  width: '400px',
  height: '200px'});
Map.add(chart);

// Outline and center the map around coordinates.

var BakersfieldLayer = ui.Map.Layer(boundary, {color: 'FF0000'}, 'Bakersfield');
Map.layers().add(BakersfieldLayer);
Map.setCenter(-118.72, 35.4, 9);

// Create a label on the map.
var label = ui.Label('Click a point on the chart to show the image for that date.');
Map.add(label);

// When the chart is clicked, update the map and label.

chart.onClick(function(xValue, yValue, seriesName) {
  if (!xValue) return;

// Show the image for the clicked date.

  var equalDate = ee.Filter.equals('system:time_start', xValue);
  var image = ee.Image(L8_TOA.filter(equalDate).first());
  var l8Layer = ui.Map.Layer(image, {gamma: 1.3, min: 0, max: 0.3,});
  Map.layers().reset([l8Layer, BakersfieldLayer]);

// Show a label with the date on the map.

  label.setValue((new Date(xValue)).toUTCString());});

  // Predefine the chart titles.

var title = {
  title: 'Bakersfield site',
  hAxis: {title: 'Date'},
  vAxis: {title: 'NDVI'},};

print(Chart.image.series(NDVI_bands, geometry, ee.Reducer.mean(), 200).setOptions(title));


texts: Dynamically changed url in my Vector layer protocol not taken into account after refresh

texts:
I have a map which I try to change dynamically its url when the user clicks a radio button.

Here is my code :

var refreshLayer = new OpenLayers.Strategy.Refresh({ force: true, active: true });
var fromProjection = new OpenLayers.Projection("EPSG:4326");
var toProjection = new OpenLayers.Projection("EPSG:3857");
var vectorProtocol = new OpenLayers.Protocol.HTTP({
        url: '/getmaindata',
        format: new OpenLayers.Format.GeoJSON({
            'internalProjection': fromProjection,
            'externalProjection': toProjection
        }),
        params: {
            "naf": "04"
        }
    });
var map = new OpenLayers.Map("map_element");
var vectorStrategies = [new OpenLayers.Strategy.Fixed(), refreshLayer];
var osmLayer = new OpenLayers.Layer.OSM();
var vectorLayer = new OpenLayers.Layer.Vector('vector_layer', {
   protocol: vectorProtocol,
   strategies: vectorStrategies
});

map.addLayers([osmLayer, vectorLayer]);

(...)

// $radio is my radiobutton
$radio.onclick = function() {
      vectorProtocol.url = "/getouterdata";
      refreshLayer.refresh({force: true});
});

When I debug with Firebug, I notice that the url is correctly changed in my layer's protocol object. However, what is called is the previous url (/getmaindata) instead of the new one (/getouterdata).
So I was asking myself if there was something I did wrong. When is the new url called ? Is it by calling the refresh() method of my Strategy.Refresh ? Or somewhere else ? Maybe the old url is cached ?


texts: How to manually geocode records or see a reference list of countries CartoDB recognizes?

texts:
I'm making a quick choropleth map that classifies countries by income level.

The in-built geocoder works mostly fine on the country names (there are also 3 letter codes in the dataset) but there are 8 or so areas that it hasn't geocoded and assigned a polygon to, for example "Isle of Man", "Andorra", "Kosovo" and "East Timor".

Is there a way for me to manually geocode these or to see what reference list of countries CartoDB is using so that I might tweak the dataset to get the geocoder to work?


texts: Show names once on a dual carriageway in Tilemill

texts:
I am trying to render a local map using OSM data with Tilemill.

The biggest issue is having road names show up twice for a dual carriageway.

I was thinking I could use oneway="1" plus another selector to have it only render on "odd" streets but could not figure a way to make this work.

Here is my css:

#osm.roads[type="trunk"],
#osm.roads[type="primary"] {
  ::label {
    text-name: [name];
    text-face-name: 'Droid Sans Regular';
    text-placement: line;
    text-size: 15;
    text-dy: 12;
    text-halo-fill: fadeout(white, 30%);
    text-halo-radius: 2.5;
    text-max-char-angle-delta: 17;
    text-min-distance: 5;
  }
}


texts: getCenterLonLat differs from init values

texts:
bounds = new OpenLayers.Bounds();
bounds.extend(new OpenLayers.LonLat(36,57));

console.log(bounds.getCenterLonLat());

I expect something  like  Object { lon=36, lat=57, ... }

I get some Object with a different value for lon.

What can I do?

Question still open.


texts: What is a good data structure and file format for describing road network geometry AND topology?

texts:
I am trying to emulate the road network from OSM or Google (or Waze, etc.) locally, so that I can run some algorithms on it.

For that, I want to represent these roads as graphs, where each edge would be a linestring, each node would be a point, and each point, etc.

Is there already a good file format for that, preferrably text-based? (JSON or CSV would be great). I have already seen some formats, such as STL, which stores geometries as indexed points, and describe the topology only by referencing these points. As such, an intersection would refer the same point from each of its incoming or outgoing roads, for example.

I am not a GIS professional, only an interested programmer, so I would prefer simplicity and light-weight-ness if possible. I don't mind devising my own custom "format" if I get to understand what type and structure of information I would be supposed to be representing in the first place.


texts: arcpy.CreateFeatureclass_management ERROR 000840 after publishing?

texts:
I have a exception in arcpy script.In my code I can get  result  of Trace  on  geometic network as group layer..next  I want  to  extract only points  and return  those  values  as FetureLayer .So  I must to  create a Featureclass and  insert the  results of trace  in to it, this  work  for me  on ArcMap but When I publish it on Arcgis Server   I  get  this  error

in CreateFeatureclass raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000840: The value is not a Workspace. ERROR 000840: The value is not a Feature Dataset. Failed to execute (CreateFeatureclass).

and  it is  my code

FC =arcpy.CreateFeatureclass_management("in_memory","M12","POINT","","DISABLED","DISABLED",spatial_reference)
 arcpy.AddField_management(FeatureClass, "TableName", "TEXT", "", "", "80", "", "NULLABLE", "NON_REQUIRED", "")
 arcpy.AddField_management(FeatureClass, "Obj_ID", "TEXT", "", "", "80", "", "NULLABLE", "NON_REQUIRED", "")
 arcpy.AddField_management(FeatureClass, "GeoJson", "TEXT", "", "", "880", "", "NULLABLE", "NON_REQUIRED", "")
 arcpy.AddField_management(FeatureClass, "LayerId", "TEXT", "", "", "80", "", "NULLABLE", "NON_REQUIRED", "")
 arcpy.AddField_management(FeatureClass, "LayerName", "TEXT", "", "", "80", "", "NULLABLE", "NON_REQUIRED", "")

Does anybody know where is the problem?


texts: ogr2ogr merge fgdb (while importing to PostGIS)

texts:
I have millions of records and hundreds of tables in several file geodatabases..

The tables name in these FGDBs are all same but just in different areas.

For instance, the Admin table, you see below, is California data (but the FGDB named as file.gdb), and there is another Admin table in New York (also the FGDB named as file.gdb) data

To have better access and use st_* functions, I was importing the FGDBs to my PostgreSQL one by one. (ultimately I am going to write a bash script to automate this)

As far as I know, once the tables are imported to PostgreSQL I can append tables using -update -append options.

So, this is my command below

ogr2ogr -update -append -progress -f "PostgreSQL" PG:"host=localhost port=5433 dbname=postgres user=postgres password=postgres" file.gdb -lco GEOMETRY_NAME=shape -t_srs EPSG:3857

but this gives me this error message

  Warning 1: Layer creation options ignored since an existing layer is
         being appended to.
  ERROR 1: ?ERROR 1: INSERT command for new feature failed.
  ?ERROR 1: Unable to write feature 1 from layer Admin.
  ERROR 1: Terminating translation prematurely after failed
  translation of layer Admin (use -skipfailures to skip errors)

Using -skipfailures seems to be dangerous because it might skip some of my data.

Will there be another solution to breakthrough this issue?


texts: OS X Yosemite: copy / paste into shell GRASS 6.4.4. within QGIS 2.14.0

texts:
I want to know How to copy/paste into GRASS 6.4.4. shell  within QGIS 2.14.0.
I'm working with OS X Yosemite and the right mouse button doesn't work neither shortcut (cmd+c/cmd+v)


texts: Looking for a SVG library for construction

texts:
Would anyone be able to suggest a good library for svg symbols for construction please? To be used in QGIS. Ideally water and pipeline network symbols.

ciao


texts: ``Qt: QUndoStack::endMacro(): no matching beginMacro()`` error on PyQGIS script

texts:
I am working on a QGIS plugin that allows redrawing a selected feature without losing its original attributes.
The plugin works properly and the process ends well but I am receiving a QT error that I can't solve: Qt: QUndoStack::endMacro(): no matching beginMacro()
This is the portion of the code:
  def run(self):
        """Run method that performs all the real work"""

        if layer is not None:

            # Check, with messages, that only one feature is selected
            if layer.selectedFeatureCount() == 0:
                QMessageBox.warning(
                None,
                'Replace Geometry Plugin',
                'There is no selection. Select one feature.')
            elif layer.selectedFeatureCount() >= 2:
                 QMessageBox.warning(None,
               'Replace Geometry Plugin',
                'The selection contains multiple features. Select only one feature.')

            else:

                # Turn off the attribute form popup
                form_config = layer.editFormConfig()
                form_config.setSuppress(QgsEditFormConfig.SuppressOn)
                layer.setEditFormConfig(form_config)

                for feature in layer.selectedFeatures():
                    global old_id
                    old_id = [feature.id()]

                def feature_added(self):

                    # Disconnect from the signal
                    layer.featureAdded.disconnect()

                    listOfIds = [feat.id() for feat in layer.getFeatures()]
                    new_feature_id = [listOfIds[0]]

                    # Calculate the new WKT
                    layer.selectByIds(new_feature_id)
                    for feature in layer.selectedFeatures():
                        new_WKT = feature.geometry().asWkt()
                        layer.removeSelection()

                    layer.startEditing()
                    layer.selectByIds(old_id)
                    for feature in layer.selectedFeatures():
                        final_id = feature.id()
                        new_geometry = QgsGeometry.fromWkt(new_WKT)
                        layer.changeGeometry(final_id, new_geometry)

                    # Remove the new feature
                    layer.removeSelection()
                    layer.selectByIds(new_feature_id)
                    for feature in layer.selectedFeatures():
                        toremove_id = feature.id()
                        layer.deleteFeature(toremove_id)
                        # Save the changes
                        layer.commitChanges()

                    # Reselect the old feature
                    layer.selectByIds(old_id)

                    # Turn on the popup if ever was on
                    form_config = layer.editFormConfig()
                    form_config.setSuppress(QgsEditFormConfig.SuppressDefault)
                    layer.setEditFormConfig(form_config)

                # Connect the layer to the signal featureAdded, so when a feature is added to the layer, the feature_added function is called
                layer.featureAdded.connect(feature_added)

                # Set the layer in edit mode
                layer.startEditing()

                #Activate the QGIS add feature tool
                iface.actionAddFeature().trigger()

Any suggestion to solve it and avoid the error?


texts:
The coastline (the data I'm using) is including the small islands (Tarakan and Laut), I wanted to create two labels for these islands. In attribute table itself doesn't have specific name, so for this, I inserted a column named 'name_hoel' in the attribute table, after setting, it shows me only the 'Laut', but NOT 'Tarakan'.

What shall I do to let both islands name appear? (at a scale 25 Mio.)

Layer Properties>Labels>Rendering and the box says: "Show all labels for this layer (including colliding labels)"---- checked! but it is still not showing.

image: 1) coastline with wanted islands (in yellow)

   2) it should look like this (please ignore its color):

   3) Layer Properties: 'placement' and 'rendering'


texts: gdal2tiles rx, ry, rxsize, rysize, wx, wy, wxsize, wysize

texts:
Reading the source code of gdal2tiles.py, one often comes across the variables rx, ry, rxsize, rysize, wx, wy, wxsize and wysize.
For example, when reading raster data, you encounter:
data = ds.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize, band_list=list(range(1, dataBandsCount + 1)))
How does rx, ry, rxsize, rysize, wxsize and wysize relate to the actual raster file being read? It's not immediately obvious from reading the documentation.
Further down from there, gdal2tiles.py does:
# TODO: fill the null value in case a tile without alpha is produced (now
# only png tiles are supported)
dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1, dataBandsCount + 1)))
dsquery.WriteRaster(wx, wy, wxsize, wysize, alpha, band_list=[tilebands])

scale_query_to_tile(dsquery, dstile, tile_job_info.tile_driver, options, tilefilename=tilefilename)

Basically the question is, what are the `r` variables and the `w` variables in relation to the source raster?


texts: Does gdal_calc only support 26 input raster files at a time?

texts:
In gdal_calc.py, input files are specified with letters A-Z.

Here are two files, input.tif and input2.tif, being averaged into result.tif:

gdal_calc.py -A input.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"

I have more than 26 files. What can I do about this?

One possibility: keep adding up files, 25 at a time, and then divide the values by the total number of input files. I think my numbers will become too large.

It looks like the source code for gdal_calc.py could be adapted to process as many images as needed, but the script is a mammoth nested looping function called doit(), so that's a bit scary. Should I really be hacking source code, or is there a better way?

My images are all single-band GeoTIFFs, by the way.

EDIT: After posting, I realized this question is somewhat of a duplicate of this rather incoherent one. Fortunately, the other question does have a good answer showing how to average more than 26 images together from the command line. The answer only works for a predetermined number of images, though. It would be nice to have a general solution that supports a theoretically unlimited number of images.


texts: How to select all points within a circle in PostGIS?

texts:
I have set of points in a PostGIS Table (Type Geometry).

Now, I need to filter points that lie in a circle with a center (C) and radius (R)!

I expect C to be of type Point or Coordinates..  R to be in Meters!

SELECT * FROM my_points
    WHERE ST_Within(point, ...)

However, I don't know how to complete it..

How can I exactly measure C and R?


texts: GDAL python scripts not importing in Linux

texts:
I'm new to both GDAL and Linux. I have installed GDAL on my Linux mint PC. I am able to use regular GDAL commands in the terminal, but when I use python based commands,such as gdal_calc.py I get the following error:

Traceback (most recent call last):
  File "/usr/bin/gdal_calc.py", line 51, in <module>
    from osgeo import gdal
  File "/usr/lib/python2.7/dist-packages/osgeo/__init__.py", line 21, in <module>
    _gdal = swig_import_helper()
  File "/usr/lib/python2.7/dist-packages/osgeo/__init__.py", line 17, in swig_import_helper
    _mod = imp.load_module('_gdal', fp, pathname, description)
ImportError: /usr/lib/x86_64-linux-gnu/libspatialite.so.7: undefined symbol: freexl_version

I get the same error when I try from osgeo import gdal in a Python session. I am presuming that this has something to do with the locations of the gdal python scripts, or something related to spatialite (although I am not sure what spatialite has to do with this command). I also have QGIS and PostGIS installed, which presumably could cause some conflicts.

I have searched numerous sources, but can't find the cause of this error. I wonder if someone might be able to point me in the right direction.


texts: GeoServer cuts parts of my raster images out and fill other ones with the black color

texts:
I have some historical maps that I would like to serve by a GeoServer.

These historical maps are GeoTIFF files that use EPSG:4326 and they have the alpha band.

I've opened them in my QGIS and they are connected correctly.

However when I serve them inside a GeoServer, with default parameters, for some reason, that I don't know, parts of my rasters are cut out and other ones are filled with the black color, in both QGIS visualization by WMS and GeoServer layer preview.

QGIS visualization by WMS:

GeoServer layer preview:

I don't understand why this event is happening, because the first and second pictures above illustrate the same GeoTIFF files, where the first picture they are showed by QGIS and the second one they are served by GeoServer using WMS.

What is wrong?

EDIT 1: Ian Turton's questions.

How did you add them to GeoServer?

I have a GeoServer 2.16.2 inside a Docker container with a volume inside my host machine where the rasters are. Then, I've created a Mosaic with default options (I did not change anything) that points to the volume with my rasters.
Is it a mosaic or individual rasters?

Mosaic. I added a new Store > ImageMosaic and I published the layer with default options.
What projection did you set?

The projection that I did set to my rasters is EPSG:4326.
What projection is the map you are viewing?

My background image is OSM and its default projection is EPSG:3857.
Does the log file have any interesting messages for you?

No, it does not, unfortunately. The layer preview of my raster inside GeoServer is equals to the QGIS visualization by WMS.

EDIT 2: Layer information inside GeoServer.

Inside this Gist I've added print-screens to illustrate the layer information inside GeoServer.


texts:
When using a reducing function such as .reduceRegion(), a provided projection is necessary.
This projection has to be provided as  CRS String, and crsTransform List. Both are optainable with :
var proj = Im.projection().getInfo()
var crs = proj.crs
var transform = proj.transform

yet the Funktion .getInfo() is not usable inside .map(), due to its correspondence with the Javascript Interface of Google Earth Engine.
To call the projection with the designated ee.Projection() functions is not possible, due to the expectations of .reduceRegions()
var proj = Im.projection()
var crs = proj.crs()
var transform = proj.transform()

how can a projection Inside the .map() function be retrieved for usage inside .reduceRegion()?
var min = Im.reduceRegion({
  reducer:ee.Reducer.min(),
  crs:crs,
  crsTransform:transform,
})


texts: Using QGIS Difference with Point and Polygon

texts:
After installing a fresh QGIS 3.10 version on a Windows computer I recently saw an unexpected behavior of the "Difference" feature from QGIS. (on one of my old computer with QGIS 3.4 it is different)
I'm trying to remove points from a layer that falls into a polygon from another layer.
After running the "Difference" model the points that falls into the polygon layer are moved to the x,y position x=0 and y=0.
Has anyone an idea why this happens. My expected output is simple a layer without the points that falls into the polygon.
I want to use the "difference" function because it should be integrated in a bigger model.

Process:


texts:
I have a .shp of lines where I'm able to calculate the sum of segments length but I'm trying to put a point at a given distance in meters.

The problem is it puts a point at every segment. When I give it the instruction to put a point at 10000m, if the segment is 9000m long, then it will put a point at the end, it will not pass to the next segment.

How can I iterate on each segments? Like:

for elem in elems():
  if (distance given = length)
    put a point
  else
    distance = distance - length
    ## then the loop continue to the next segment with a new distance

Or maybe another way to do this?


texts: Cannot Compile and add GDAL library for Visual C++ project in Visual Studio IDE

texts:
I wanted o create a project using the GDAL library and its API for C code.
I can create an executable for this code by entering the following on the command line in a Linux Machine:
gcc -o myprogram myprogram.c -lgdal -lm, but I also wanted to create an executable to run the same program on windows. I found the steps found in the following post: https://lazyadminblog.blogspot.com/2018/07/how-to-build-gdal-from-source-for.html
Steps:

Download from http://www.gdal.org the source package for GDAL
Extract the .zip file and generate a _vcxproj.bat script in Powershell with parameters by typing .\generage_vcxproj.bat 14.0 64 gdal_vs2015
Execute the file created "gdal_vs2015.vcxproj"
Build in Visual Studio by typing CTRL + SHIFT + B
Create a VS Code Studio Project and open the properties and add the following directories in the "Include Directories" and "Library Directories" the following paths: $C:\Users\user_name\gdal-3.0.4 $C:\Users\user_name\gdal-3.0.4\gcore $C:\Users\user_name\gdal-3.0.4\apps $C:\Users\user_name\gdal-3.0.4\port $C:\Users\user_name\gdal-3.0.4\ogr
Add in Linker -> Input -> Additional Dependencies the following path: "C:\Users\user_name\gdal-3.0.4\gdal_i.lib"
Create executable of the C Code

The process ended creating a .exe file, which I understand it as being the executable of the program for a Windows Platform. However when I execute the program I got an Error along this lines: System error: The program can't start because MSVCP140D.DLL is missing from your computer. Try reinstalling the program to fix this problem. I did reinstall gdal, but it did not solved the problem. Did anyone confronted this problem, and found a solution to it?


texts: SVG markers made in InkScape appear transparent in QGIS

texts:
Using QGIS 2.14.3 on Ubuntu 16.04, I created a custom icon using inkscape for QGIS. You can find it here (Update: that icon is now fixed). In inkscape, in my filesystem, it displays as:

However in QGIS, it displays as

The colors are transparent. How can I fix this?


texts: How to install shapely for python 2.6 (Mac) NOT 2.7?

texts:
I'm trying to use Polygonizer plugin inside QGIS, but it needs Shapely to work. I've installed Shapely (pip install Shapely) but that's for python 2.7. QGIS works with python 2.6 (on Snow Leopard) and I can't find previous installers for it. Maybe I'm missing something?

Thanks.


texts: Osm2pgsql, permission denied when creating node cache file

texts:
This is probably a dumb issue but I'm kind of stuck while trying to import a big data set into an existing database with osm2pgsql.

More specifically, I'm getting a permission denied error when using the --flat-nodes flag followed by /mnt/nodes/flatnodes.bin, whereas typing just flatnodes.bin does not throw this error. I'm using four disks, one is for the system disk, which I don't want to use neither for the database, the flatnodes.bin nor the downloaded data set due to limited size (50GB); the bigger disk is used to store the PostgreSQL database and it is here that I intend to import the map into; a third one, which is used to store the flatnodes.bin files used during the import process; and finally a fourth disk, which is where the downloaded data (planet-latest.osm.pbf) is stored. These last two are deleted at the end of the import process. Here's what I'm typing into the command line:

postgres@map-server / $ osm2pgsql -c -d postgres --slim -C 30000 --flat-nodes /mnt/nodes/flatnodes.bin /mnt/data/planet-latest.osm.pbf

Here's what I'm getting in return:

osm2pgsql SVN version 0.88.1 (64bit id space)

Using built-in tag processing pipeline
Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
Setting up table: planet_osm_line
Setting up table: planet_osm_polygon
Setting up table: planet_osm_roads
Allocating memory for dense node cache
Allocating dense node cache in one big chunk
Allocating memory for sparse node cache
Sharing dense sparse
Node-cache: cache=30000MB, maxblocks=3840000*8192, allocation method=11
Mid: loading persistent node cache from flatnodes.bin
Failed to create node cache file: Permission denied
Error occurred, cleaning up

If I do it as root, I get this in return:

Osm2pgsql failed due to ERROR: Error: Connection to database failed: FATAL:  role "root" does not exist

Doing all this on Ubuntu Bionic with PostgreSQL 10 and PostGIS 2.5.0 installed.


texts:
Kindly advise or assist, I want to achieve following in overpass turbo:

Map/display transmission lines (power=line) and distribution lines (power=minor_line) of a specific area in map. I used below query but it doesn't correctly work (it only displays circles in map, not powerlines):

area[name="Bonn"]->.a;

(node[building](area.a);
  way[power=cable](area.a);
  rel[building](area.a);
);
out center;

My query should display powerlines in map as below figure.


texts: Iterating through records specific number at a time using ArcPy?

texts:
I'm working with a table that has over 180,000 records.  My code searches through the table and pulls field values (meeting a SQL statement) into a list.  This list is used as part of a Search by Attribute SQL statement on a feature class using the IN parameter.  For example, SELECT from fc WHERE field IN('value1', 'value2', etc).

Only posting the relevant code here.  The problem is that the initial query might return as many as 180k+ records and thus the resulting where clause is a massive list.  I've determined that if the list of values (GUIDS as strings) gets over 20,000 or so, the script fails with an "Invalid Expression" error.

The only thing I can think to do is modify this code so that it runs the first query on the table no more than 10,000 rows at a time.  That should keep the IN clause manageable.  Once that list is used to select features with matching field values in the fc, the next 10k would run as "Add_to_selection."

However, I'm having a hell of a time figuring the best way to manipulate this bit of code to do it that way.

Summary:  Both the table and the fc have a field called "Main_ID".  The code searches for those values in the table based on a user query, puts them into a list and searches on the "Main_ID" field of fcs in the mxd using a sql expression with an IN clause.  I think this line needs to be altered to just do this maybe 10k records at a time:

set([row[0] for row in arcpy.da.SearchCursor(linkTable, linkKeyField, sqlExp)])

Any help?

linkKeyField = "Main_ID"

if len(layers) > 0:
    # Set the SearchCursor to look through the selection of the linkTable
    sourceIDs = set([row[0] for row in arcpy.da.SearchCursor(linkTable, linkKeyField, sqlExp)])
# Add DBMS-specific field delimiters
fieldDelimited = arcpy.AddFieldDelimiters(arcpy.Describe(linkTable).path, linkKeyField)
# Add single-quotes for string field values
valueList = ["'%s'" % value for value in sourceIDs]
# Format WHERE clause in the form of an IN statement
whereClause = "%s IN(%s)" % (fieldDelimited, ', '.join(map(str, valueList)))
arcpy.AddMessage("SQL Clause: {}".format(whereClause))
for lyr in layers:
    if len(arcpy.ListFields(lyr, "Main_ID")) > 0:
        # Process: Select Layer By Attribute
        arcpy.AddMessage("Querying related records in {0}".format(lyr))
        arcpy.SelectLayerByAttribute_management(lyr, selectType, whereClause)


texts: Re-projecting shapefile from an unknown projection

texts:
I have been given a shapefile which is supposedly in a NAD1927_Zone16 projection. I need to align it with a raster image in a Lambert Conformal Conic (NAD 1927) projection.

I've tried to use the Project and Define Projection tools to verify that this shapefile is in NAD1927_Zone16, but when I compare to other data it is no where near where it should be (leading me to believe that the projection indicated on the shapefile is not correct).

In addition, I've edited the Standard Paralell and Latitude of Origin and Central Meridian values as per:
https://geonet.esri.com/thread/9764

And I've also tried loading in a custom projection from a .prj file downloaded from the following site:
http://spatialreference.org/ref/esri/102002/

I haven't done much with projections or transformations before, but I really didn't expect it to be this complicated when I agreed to the task.


texts: Pivot a polyline at vertex or distance using ArcObjects C#

texts:
I hope you can help me with that one. I'm using ArcObjects with C# trying to rotate only one part of a polyline at a vertex or distance along the polyline. The polyline is zAware and has values at regular spacing.

The only successful way I've found so far is splitting the polyline at the vertex first, then using ITransform2D.Rotate to rotate one of the resulting features by a specified angle, and finally merging both features again afterwards.

All this takes a long time and I had many lines to code. Unfortunately I haven't been able to find a quicker solution. But if any of you know about a more clever way, please tell me!

Thank you very much!


texts:
I am a student of geophysics and I am studying geographic information systems (ArcGIS) and we have been given a project. As I do not have much information about GIS capabilities can you please suggest me a topic for my project which is related to geology or geophysics?

For example, my fellow classmates are doing project on they creating a tectonic map which will also show where earthquake occurs.


texts: PostgreSQL function with BOX3D error "doesn't start with BOX3D"

texts:
I created function search nearest point (pois) by using pg_routing as following:

CREATE OR REPLACE function find_pois(inn_point text, distance int) RETURNS float[] AS  --point(x y)
$$
DECLARE
    source_id INTEGER;
    target_id INTEGER;
    row_data RECORD;
    i integer;
    retVal float[];
    temp_length float;
BEGIN
 -- search edge of source
 SELECT into source_id source ,
      ST_Distance(the_geom, ST_GeometryFromText(
       'POINT(587314.857692619 2325407.8399756)', -1)) AS dist
    FROM ways
    WHERE the_geom && setsrid(
       'BOX3D(587214.857692519 2325407.8399756, 587414.857692619 2325507.8399756)'::box3d, -1)
    ORDER BY dist LIMIT 1;

 -- get all point in 500m radius
 FOR row_data IN SELECT b.gid, x(the_geom), y(the_geom) FROM station b
   WHERE ST_DWithin(the_geom, 'POINT(587314.857692619 2325407.8399756)', 500.0)
    LOOP
         -- with each point, search shortest part
         -- 1. find taget_id

         SELECT into target_id target ,
          ST_Distance(the_geom, ST_GeometryFromText(
           '''POINT(' || row_data.x || ' ' || row_data.y || ')''', -1)) AS dist
        FROM ways
        WHERE the_geom && setsrid(
           '''BOX3D(' || row_data.x-100 || ' ' || row_data.y-100 || ', ' || row_data.x+100 || ' ' || row_data.y+100 || ')'''::box3d, -1)  --cast((cast(row_data.x as float))-100 as text)
        ORDER BY dist LIMIT 1;
    -- 2. tim duong di ngan nhat tu diem nguon toi diem dich
    SELECT --rt.gid, ST_AsGeoJSON(rt.the_geom) AS geojson,
        into temp_length sum(length(rt.the_geom))
        --ways.gid
        FROM ways,
           (SELECT gid, the_geom FROM dijkstra_sp_delta( 'ways', source_id, target_id, 100) ) as rt
            WHERE ways.gid=rt.gid;
       retVal := temp_length;
    END LOOP;

 return retVal;

END;
$$ LANGUAGE plpgsql;

Everytime i try excute this function, i have error:

ERROR:  BOX3D parser - doesnt start with BOX3D(
LINE 1: ... || ' ' ||  $4  || ', ' ||  $5  || ' ' ||  $6  || ')'''::box...
                                                             ^
QUERY:  SELECT target , ST_Distance(the_geom, ST_GeometryFromText( '''POINT(' ||  $1  || ' ' ||  $2  || ')''', -1)) AS dist FROM ways WHERE the_geom && setsrid( '''BOX3D(' ||  $3  || ' ' ||  $4  || ', ' ||  $5  || ' ' ||  $6  || ')'''::box3d, -1) ORDER BY dist LIMIT 1
CONTEXT:  PL/pgSQL function "find_pois" line 27 at SQL statement

I guest there are some error on string in BOX3d but cannot fix
Please help me.


texts:
We are carrying out several research projects involving volume calculations using UAV imagery for which we need survey grade GPS / GNSS accuracy for our Ground Control Points. Right now we are working on very limited budget and are beg/borrow/steal-ing access to Dual channel commercial equipment from government and industry partners. We are quickly hitting a point where we need to operate our own equipment. We do not need RTK, and can post process our corrections.

I have read about Reach and Piksi, but I do not need to go that low cost, and also want something fairly turnkey. I am not looking to integrate any of this into the UAV's, just capture accurate readings for my GCPs.

Is there a solution I am not turning up that is affordable ($1500-$3000) that does this?
If not does anyone have experience with Reach or Piksi in the field?
Has anyone evaluated the above against other high end field equipment?

We are working in Arctic and high Arctic conditions so some robustness and reliability is required.


texts: "Impossible to engage a segment in the current layer"

texts:
How can I fix this ? My problem is the red "x".

this is the screen http://imageshack.com/a/img540/8232/o6zUnB.png


texts: Applying different equation based on conditional statement for raster value in R?

texts:
I am trying to find the best way to apply a set of equations to a raster conditional on values in a separate raster. For example, suppose r1 has the values of interest and raster 2 (r2) is a raster with three values, 1, 2, and 3.

r1 <- raster(nrow=5, ncol=5)
r1[] <- runif(ncell(r1)) * 10
r2 = raster(nrow=5,ncol=5)
r2[] = round(runif(ncell(r2),1,3))
r3 = 0/r2 ## holder for values

I have modified code  that applies a conditional function and applies it to a matrix, then assigns the values of the matrix to the raster.

Con=function(condition, trueValue, falseValue){
 return(condition * trueValue + (!condition)*falseValue)}

fun1 = function(x){x*1}
fun2 = function(x){x*2}
fun3 = function(x){x*3}
m = as.matrix(Con(r2 == 1,fun1(r1),Con(r2 == 2, fun2(r1), Con(r2 == 3,fun3(r1),0))))
values(r3) = m

This works, but (a) the conditional statement is highly unreadable, and (b) I am concerned that is will be prone to error with larger rasters.

Does anyone have a better way of doing this?


texts:
I try to follow the step in How to import qgis.core to Python app? to import qgis.core with Python Shell. But when I type

set PYTHONPATH=C:\Program Files (x86)\Quantum GIS Wroclaw\apps\qgis\bin

it report

SyntaxError: invalid syntax

and highlight PYTHONPATH


texts: 3D model with added annotations for phone app

texts:
This might be too broad of a question; I have some 3D models that I am viewing in CloudCompare. I can add a label to a desired point by using shift + right click and then change the title of the point to say whatever I want. I then am wanting to export the model with the added annotations and open it up in a phone app. I have been using EMB3D, but am open for suggestions (would like to stay in the free realm). I have not been able to open the file into EMB3D because the file types are not supported.

I am wondering if anyone has any suggestions or resources that I could read through for something like this. Basically what I want is to simply add annotations to desired locations on my 3D model and view it on a phone.


texts: Adding more than 5000 paths string in ElevationProfile using ArcGIS API for JavaScript?

texts:
I'm creating a simple ElevationProfile using the ArcGIS API for JavaScript.

In the URL it's described as

MaximumRecords: 5000

Its allows only 5000 polyline array paths but I have more than 5000 polyline array paths (ex : 12,000).

Is it possible to send the more than 5000 polyline path? If not, can I Compress the large polyline array into below 5000?

Give suggestion for the above and how to proceed?


texts: Is it possible to render TIFF image with QGIS Server?

texts:
I would like to know if QGIS Server is able to render a TIFF Image ?

Actually, I have a QGIS Server running, this server serve correctly my QGIS project to an OpenLayers (java script library) frontend. Images are requested by OpenLayers in PNG, and every things works fine.

Now I would like to get some specific layers in TIFF image containing float32 value. These layers are DEM layers, and I need to get the precise altitude on each pixel (and not a colorscale rgb value). I can't find anywhere in QGIS an option to force this layer to be served in TIFF Float 32 (which is the raw type of my layer's image).

When I look to the GetCapabilities of QGIS Server, I can read :

<GetMap>
<Format>image/jpeg</Format>
<Format>image/png</Format>
<Format>image/png; mode=16bit</Format>
<Format>image/png; mode=8bit</Format>
<Format>image/png; mode=1bit</Format>
<Format>application/dxf</Format>
<DCPType>...</DCPType>
</GetMap>

So is it possible to add a format in this list ?