texts: Building SQL statement in javascript and casting variables as strings

texts:
I'm working with Postgres 8.3/PostGIS which is storing some polygons. I'm writing a javascript function which will take in a parameter of a watershed id from another function. The javascript will build up a sql statement and then send that along to a jsp script which will send that to the database and return some geojson which will then be output into Openlayers. The problem is that the geojson is not being returned because I think I have a problem with my select statement and some type (integer/text) problem.
The table holding the polygon geometry contains a column set as varchar(50) but it contains some id's which have only numbers, and some id's which contain numbers and string characters. The function attempts to return a geojson from the database of a polygon of a given id. I've left off the last end of the function because I think my problem is somewhere in the building of the sql statement.

function displaySelectedShed(shedId) {
    var json_url = "04-sql-to-json.jsp?sql=";
    var sql_geom = "st_transform(ST_Force_2D(the_geom),900913)";
    var sql = "select st_asgeojson(" + sql_geom + "), id ";
        sql += "from medford.subshed2 ";
        sql += "where subshed2.id = " + shedId +"";
        json_url += encodeURIComponent(sql);
 }

The select statement above looks like this:

select st_asgeojson(st_transform(ST_Force_2D(the_geom),900913)), id from medford.subshed2 where subshed2.id = 01616500

and I can test this by putting this:

http://localhost:8080/spatialdbtips/04-sql-to-json.jsp?sql=select%20st_asgeojson%28st_transform%28ST_Force_2D%28the_geom%29%2C900913%29%29%2C%20id%20from%20medford.subshed2%20where%20subshed2.id%20%3D%2001616500

into a browser window. But the error I get is this:

Error: operator does not exist: character varying = integer

I know this has probably something to do with the fact that my shed id is a number when inserted as a variable into the sql statement, but I've tried all sorts of cast combinations to get it to be a string, since the values in the database are stored as such, but I can't get them to work. I've tried both the CAST as well as the  :: such as

sql += "where subshed2.id = " + shedId +"::text";

but I can't get it to work. Does anyone have any suggestions to get my watershed id numbers (the shedId variable) to act like string in the sql where clause, and how to build it up in a way I've shown in javascript?


texts:
I have shapefiles containing coordinate system WGS84, when i try to upload google satellite from Open Layers Plugin, it automatically changes its coordinate system to WGS84/world mercator EPS3395 (the  projected one) but fits exactly in the same position, but when i try to change its coordinated system its position got miss matched, the same thing happens when i try to do compose a map or/and export it as data

I am really new to QGIS.


texts: Make QGIS python plugin for both versions 2.x and 3.x?

texts:
I am in the process of migrating a QGIS python plugin from QGIS 2 to QGIS 3, and browsing various resources.

It's not clear if it's possible to have the plugin compatible with both versions or if it's necessary two handle to plugin versions.

The problem I hit so far is how to manage PyQt import (PyQt4/PyQt5) ?


texts: How can rewrite this query for using with commonly used lat-long values

texts:
I have this query to get the nearest road vertices from a bus stop point. Ho can a re-write this query to pass commonly used lat-long values like in google maps instead of  715127 957625

SELECT * FROM tvm_city_busroute_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText('POINT(715127 957625)',4326)
LIMIT 1;


texts:
I have two geometries and I want to write the information from one geometry into the other geometry, but only where they intersect. Sometimes, a geometry is intersected with to different geometries (see my screenshot). The yellow geometry is intersected by two geometries (thick black line). The resulting table should contain the yellow geometry twice, with the two different information from the intersecting geometries.

I tried

CREATE TABLE myTable AS
SELECT g.a, g.b, g.geom, t.c, t.d
FROM green_geometries g
LEFT OUTER JOIN thickBlackLine t
ON ST_Intersects(ST_PointOnSurface(g.geom), t.geom);

However, the resulting table only contains the yellow geometry once, with only one information from the intersecting geometry added to it.

*edit: thickBlackLine t is a a table, containing many polygons. In the screenshot here you only see the outline (thickBlackline) of the polygon.
Maybe on the screenshot here it becomes more clear. I am trying to perform a Spatial Join on polygons, not line and polygon(s)!


texts: Why CQL DWITHIN returns different result than PostGIS DWITHIN?

texts:
I have 2 tables in PostGIS that I publish to geoserver. both have same SRS (2039).

When I do the following query in PostGIS , I get 69 results :

 SELECT a.gid,b.gid
 FROM Table_A a INNER JOIN Table_B b
 ON ST_DWithin(a.geom, b.geom,1);

But when I do the same in CQL , I get error that result contains more than 1000 features.

This is my CQL :

http://localhost:8080/geoserver/wms?
LAYERS=a,b&STYLES=&FORMAT=image/png&SERVICE=WMS&VERSION=1.1.1&
REQUEST=GetMap&EXCEPTIONS=application/vnd.ogc.se_inimage&
SRS=EPSG:2039&
BBOX=589042.42768447,4914010.3926913,609134.78143081,4928059.3431623&
WIDTH=512&HEIGHT=358&
CQL_FILTER=
DWITHIN(geom,collectGeometries(queryCollection('b','geom','INCLUDE'),1,meters);
EXCLUDE

The Exclude at the end is because I don't want to see layer b as result.

Any explanations on this ?

* UPDATE *

Layer A and Layer B are published to Geoserver , The dimensions in the request I get by default when I use PREVIEW LAYER (Layer A are lines and B are polygons)


texts: How to effectively highlight small polygon areas in small scale maps

texts:
I've recently been working on a thematic map which shows a large area at small scale, containing a number of different polygon layers showing different feature classes (unfortunately I can't share the draft map). One thing I'm struggling with on this map is how to effectively highlight the presence of small area polygons. These tend to get lost among the larger area polygons, but they are important and need to be immediately noticeable on the final map.

I've tried a number of techniques for highlighting these, but without success:

Using multiple inset maps at a larger scale isn't an option due to the number of these small area polygons on the map. The final product also needs to fit on a single page, so space is a premium.
Labeling them also isn't an option, again due to the number of these problematic features.
I've tried adding a larger 'halo' style border under these polygons, as a way of making them visually larger while still attempting to give a good indication of their actual size. It's probably the closest I've come, but I'm not that happy with how this looks (it's more cluttered and confusing then I'd like).
I also tried swapping the feature symbol from a polygon to a point symbol when the feature is below a certain size threshold - but I think this makes for a less readable map at first viewing. I'd also like to have some way to give an indication of the actual size of the smaller features still.

So now I'm after tips or tricks in how to visually increase the prominence of these features so that they stand out more on the finished map.


texts: Fill QLineEdit with coordinates by click on QGIS canvas

texts:
With a plugin I can click a button and activate a QgsMapTool, defined as:

class PointTool(QgsMapTool):

def __init__(self, canvas):
    QgsMapTool.__init__(self, canvas)
    self.canvas = canvas

def canvasReleaseEvent(self, event):
    #Get the click
    x = event.pos().x()
    y = event.pos().y()

    point = self.canvas.getCoordinateTransform().toMapCoordinates(x, y)
    coordinates = str(point).replace('(','').replace(')','')
    coord_x = float(coordinates.split(',')[0])
    coord_y = float(coordinates.split(',')[1])

    QgsMessageLog.logMessage(str('{0:.1f}'.format(coord_x)) + ' ' + str('{0:.1f}'.format(coord_y)),'coords')

It shows me the coordinates in a tab 'coords' of the Log Messages Panel.

The main function I use is within the class class mask_n_functions(QDialog, Ui_InputMask):. And there I call the QgsMapTool as follows:

self.xy.clicked.connect(self.setToolActive)

def setToolActive(self):
    tool = PointTool(iface.mapCanvas())
    iface.mapCanvas().setMapTool(tool)

Now I want to fill two QLineEdit in class mask_n_functions(QDialog, Ui_InputMask): with coord_x and coord_y of class PointTool(QgsMapTool)).

Where do I have to call self.lblx.setText(coord_x)?

Besides, the coordinates from canvas are rounded:

x: 620832.96 -> 620833.0

y: 5713645.90 -> 5713650.0 !!! Why?


texts: Create a filter for an attribute having multiple values

texts:
I am new to WFS and am trying to create a filter to get all the features that have their basename = "XYZ". The structure of a single record is as:

<gml:featureMember>
        <core:AdminPlace gml:id="ADMN_21493454">
            <core:featureNames>
                <core:FeatureName>
                    <core:basename>ABC</core:basename>
                    <core:languageCode>ENG</core:languageCode>
                </core:FeatureName>
                <core:FeatureName>
                    <core:basename>DEF</core:basename>
                    <core:languageCode>ENG</core:languageCode>
                </core:FeatureName>
                <core:FeatureName>
                    <core:basename>XYZ</core:basename>
                    <core:languageCode>ENG</core:languageCode>
                </core:FeatureName>
            </core:featureNames>
            <core:adminLevel>4</core:adminLevel>
            <core:dummyAdmin>N</core:dummyAdmin>
         </core:AdminPlace>
    </gml:featureMember>

I am trying to use the XML input for POST to WFS as:

<?xml version="1.0" encoding="UTF-8"?>
<wfs:GetFeature service="WFS" version="1.0.0"
  xmlns:tns="http://www.abc.com/core"
  xmlns:wfs="http://www.opengis.net/wfs"
  outputFormat="GML3.1">
   <wfs:Query typeName="tns:AdminPlace">
    <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml">
       <ogc:PropertyIsEqualTo>
         <ogc:PropertyName>tns:featureNames/tns:FeatureName/tns:basename</ogc:PropertyName>
        <ogc:Literal>Ellerth</ogc:Literal>
        </ogc:PropertyIsEqualTo>
    </ogc:Filter>
  </wfs:Query>
</wfs:GetFeature>

but I am not able to get the result. The error message thrown is:

There was an error processing the event:GetFeatureException: Error while retrieving features: could not resolve property: featureNames.featureName.


texts: What is the method for mapping 200k records on Google Maps API via MYSQL db?

texts:
I have a simple Google Maps api linked to MYSQL. I am trying to map a mysql table of 200k records. I have a default display of this data set as heatmap, so am not trying to map 200k markers. Even still, the app freezes, but I am unable to get a specific error.

Is there an another method I should apply? (i.e., data compression)


texts: WMS with EPSG:3004 not returning tiles in LeafletJS

texts:
I am consuming an Italian WMS in a CRS with an EPSG code of 3004 in QGIS. It works somehow (Needed to change the CRS from 4326 to 3004), yet the definition in the GetCapabilities document is strange (states to be 4326 but uses coordinates in projected CRS).

I am now trying to create a webmap with the same WMS using leaflet and Proj4leaflet. All the network traffic towards the WMS returns 500 as status codes:

Example URL:

http://maps1.insula.it/common/cache/ramses/?service=WMS&request=GetMap&layers=ramses&styles=&format=image%2Fpng&transparent=false&version=1.1.1&width=256&height=256&srs=EPSG%3A3004&bbox=2312320.0789986267,5034917.271996673,2312448.078998631,5035045.271996661

If I replace the BBOX coordinates with a BBOX from QGIS it returns an image:

http://maps1.insula.it/common/cache/ramses/?service=WMS&request=GetMap&layers=ramses&styles=&format=image%2Fpng&transparent=false&version=1.1.1&width=256&height=256&srs=EPSG%3A3004&bbox=2310927.24273836985230446,5034502.11386785749346018,2311008.92086785659193993,5034583.79199734423309565

The bbox of the WMS is defined as:

2288384.0789999999105930,5011877.2719999998807907 : 2336386.6389999999664724,5057879.8320000004023314

Here is the non-working codepen. The releant part is this definition:

    var crs = new L.Proj.CRS("EPSG:3004", "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +no_defs",
    {
        resolutions: [
            8192, 4096, 2048, 1024, 512, 256, 128,
            64, 32, 16, 8, 4, 2, 1, 0.5
        ]
    }),
    map = new L.Map('map', {
        crs: crs,
    });

  L.tileLayer.wms('http://maps1.insula.it/common/cache/ramses/', {
    layers: 'ramses',
    format: 'image/png',
    maxZoom: 14,
  minZoom: 0,
}).addTo(map);


texts:
I'm trying to set symbology on a point layer in QGIS using a font marker.  I would like to use a font that contains symbols rather than normal characters.

The dropdown picker for font family seems to display the text using the style of the font.

This could be useful when trying to write a Word document or something.  However, it becomes really confusing when trying to pick amongst several symbolic fonts.

Is there a setting in QGIS that will display the font names using just a plain font without using the font itself to display it? I'm looking for the ability to read the fonts within the dropdown.


texts: How to batch convert state plane coordinates to lat long with multiple SPC Zones?

texts:
I'm somewhat of a GIS novice and am trying to see if this is possible to save me some time with further investigating. I have about 300 records in a CSV file with the fields SPC Zone, ECOORD, NCOORD, CRDUNITS, and a point identifier.
I was wondering if there is any possible way to batch convert these points (all with differing SPC Zones) to lat long in QGIS or some other application?


texts: What's difference between Topology and Check Validity in QGIS?

texts:
I would like to know what the difference is between the Topology Checker and Check Validity tools?

Is there anyone who knows of a good document or tutorial that explains functions of these tools? One that, if possible, has figures of vector topological errors.


texts: How to include streaming media into an OpenGeo system to monitor live-tracked animals?

texts:
I am wondering how to include streaming media (video), or camera feeds, into an OpenGeo system to monitor live-tracked animals.

I am currently doing so with different applications, but I want to do it in a whole system (open and free), preferably under OpenGeo suite and QGIS, or other open possibilities (may VLC be an solution?). For media and GPS tracking, we are currently using cheap Chinese "dashboard cams" and pet GPS collars, since we can catch these animals almost daily (that's ethology study about domesticated badgers), but we have to postprocess all the info and manually add data in the whole database, when it could be done automatically.

I think we have 2 possibilities, being the 2nd the preferred:

Streaming and server-side recording, which is safer since data loss is minor;
On-collar recording, with on-demand live camera feeds plus server-side recording (somehow as geotagged shots/stills, an option if we do m-jpeg feeds, but not the best, I think).

I have read some stuff, since I have very limited knowledge on these topics. Open streaming servers are a first step, but I am not just wondering if it is the best match, but if it is the right schema I am thinking about.

As example, one really nice, propietary option is Skyline that merges everything in a whole platform (but for UAV military purposes).

Ideas?


texts:
I have some data stored in a CSV file containing name and date ranges.

For example:

ID | Name | Start | End   |  X  |  Y
-------------------------------
1  | Test | 2014  | 2015` | ... | ... |
2  | Test | 2014  | 2016` | ... | ... |
3  | Test | 2020  | 2025` | ... | ... |

I join this data based on the ID column with shapefiles and thus get the information coupled with the shape data.

I can assign different colours and write definition queries to make the data stand out and filter the date range but I would like to split the data over multiple layers in order to show shape data for year 2014,2015,... etc. in ArcReader.

I have a programming background but no prior experience of using ArcGIS Desktop.

How can I do this?


texts: Line segment detection QGis orfeo toolbox

texts:
I would like to extract linear features from a high resolution DEM.

I have found that a "line segment detection" tool exists in the Orfeo toolbox in QGis.

I could not find easy to understand, comprehensive documentation with examples on how to use this tool.

I tried to run the tool, initially I got an error message which referred to an invalid operand type.

Then I downloaded an OTB folder for windows, tried to re-run the tool but now I get this error message: "OTB folder is not configured. Please configure it before running OTB algorithms".

Anyone can help?


texts: Terminating Python Script in QGIS console?

texts:
How do I terminate any script in QGIS python console window?

ctrl + c does not work


texts: Splitting lines table using point table in PostGIS

texts:
My problem is that I want to get the shortest route between multiple point pairs which are stored in a table (PostGIS).

In PostGIS I also have a network created using pgr_createTopology on a table containing about a 1000 MultiLineString rows. It's my understanding that to calculate routes between points with pg_routing, you need to either accept that the route goes to and from the nearest network node or you need to work around that limitation by adding new nodes to the network which corresponds to the stops you want to calculate on.

To add new nodes to my network I've tried tried to follow the answer from @alpha-beta-soup in this tread
 - Split lines into non-overlapping subsets based on points
Before following the answer I used ArcMap to snap the points to the lines. I then imported the two shapefiles (network and radars both with geometry columns "geom") into PostGIS using the PostGIS Shapefile loader (2.1).

When using alpha_beta_soup's code I get a table with all the right columns, but there are zero rows in it, I've tried running only the first part of the code, but I still only get a table with no content. This is the first part of the code:

CREATE TABLE tmp_lineswithstops AS (
WITH subq AS (
    SELECT
    ST_Line_Locate_Point(
        roads.geom,
        ST_ClosestPoint(roads.geom, radars.geom)
    ) AS LR,
    rank() OVER (
        PARTITION BY roads.gid
        ORDER BY ST_Line_Locate_Point(
            roads.geom,
            ST_ClosestPoint(roads.geom, radars.geom)
        )
    ) AS LRRank,
    ST_ClosestPoint(roads.geom, radars.geom),
    roads.*
    FROM network AS roads
    LEFT OUTER JOIN radars
    ON ST_Distance(roads.geom, radars.geom) < 0.0001
    WHERE ST_Equals(ST_StartPoint(roads.geom), radars.geom) IS false
    AND ST_Equals(ST_EndPoint(roads.geom), radars.geom) IS false
    ORDER BY gid, LRRank
)
SELECT LR FROM subq
);

The two tables look like this:

Does anyone have a clue as to what i'm doing wrong?


texts: Handling, comparing, passing datetime variables with ArcPy

texts:
There are droves of questions on here about handling dates/times, and I've read through volumes of them, but still don't understand what I am doing wrong.

I am trying to run a simple Python script in ArcGIS Pro but keep encountering errors... mostly having to do with data type conflicts. What I want to do is look through a field in a feature class table to find the most-recent (max) date, and then compare that with another variable date.

Here is what I have:

import arcpy
import datetime
import time

arcpy.env.workspace = "C:\Data\Reporting.gdb"

ptFC = "C:\Data\Reporting.gdb\RMI_pts"
listDates = []

cursor = arcpy.da.SearchCursor(ptFC, "ReportQueryDate")

for row in cursor:
    listDates.append(row)

maxDateTuple = max(listDates)
maxDateDT = datetime.datetime(maxDateTuple[0:6])

if maxDateDT < datetime.datetime.today():
    print (maxDateTuple)
else:
    print ("doh!")

The code above returns an error on line 18 stating that in integer is required but it got a tuple. The actual error is:

  Traceback (most recent call last):
    File "", line 16, in
  TypeError: an integer is required (got type tuple)

I don't know if my approach is correct: What is the best way to return a 'max' date from the table and then compare it to today's date?


texts: Count number of points outside of buffer

texts:
I'm trying to count the total points outside of a buffer zone using PostgreSQL. The intention is to find all of the schools that do not fall within a 2-mile radius of libraries.
I've either gotten results that return the total number of schools or nothing at all.
What is wrong with my query?
SELECT ST_Buffer(mdlibrary.geom, 3218, 'endcap = round join = round') AS librarybuffer, mdlibrary.anchorname AS libraryname, count(mdschool.*) AS totalschools
FROM mdlibrary, mdschool
WHERE mdschool.geom IS NULL
GROUP BY mdlibrary.geom, mdlibrary.anchorname


texts: how to install QGIS 2.2.0 on mac OS 10.6.8?

texts:
Does someone know how to install the latest QGIS on a mac os 10.6.8? The dmg package available on http://www.kyngchaos.com/software/qgis does not run on my OS version.
I was only able to find a developement version of QGIS from Sept 2013 which is not fully stable.
Thanks


texts: Coalesce field "X" if field "Y" is a duplicate

texts:
(QGIS 3.4) The image below will help explain my goal:

I'm working out how to coalesce the attributes of the field "ID" if the field "X, Y" is a duplicate.

I have highlighted some duplicate coordinates in red. My expression will create a new text field and chain all attributes in the "ID" field, delimited by a pipe (|).

For example, where I have labelled features with blue text 1,2,3, the new field will write

Special Protection Areas | Scheudled Monuments | SSSI | AONB | Special Areas of Conservtion | Scheduled monuments | Scheduled Monuments | SSSI | AONB

I will omit duplicate values later in excel.

Here's what I have so far...if( "X, Y" = "X, Y" , coalesce(  "ID" ),' ') but this just duplicates the "ID" field without any coalescing.


texts:
I am following this example

https://docs.geoserver.org/latest/en/developer/programming-guide/wps-services/implementing.html

To create a custom GeoServer WPS. I am trying to build the project using

mvn clean install

But the build fails with the following error:

[ERROR] Failed to execute goal on project hello_wps:
Could not resolve dependencies for project org.geoserver:hello_wps:jar:2.8-SNAPSHOT:
Failed to collect dependencies at org.geotools:gt-process:jar:14-SNAPSHOT:
Failed to read artifact descriptor for org.geotools:gt-process:jar:14-SNAPSHOT:
Could not transfer artifact org.geotools:gt-process:pom:14-SNAPSHOT from/to boundless (https://repo.boundlessgeo.com/main):
Transfer failed for https://repo.boundlessgeo.com/main/org/geotools/gt-process/14-SNAPSHOT/gt-process-14-SNAPSHOT.pom:
Unknown host repo.boundlessgeo.com

I have tried to change the GeoTools version to the latest version 24-SNAPSHOT, but the same error occurred.


texts: How to drop ArcSDE connection using Python scripting?

texts:
I want to be able to control connections to ArcSDE via python. Here is my script so far using ESRI's example.

import arcpy
from arcpy import env
env.workspace = r"SDEConnection"
fcList = arcpy.ListFeatureClasses()
print str(fcList) + "\n"
env.workspace = ""
arcpy.ClearWorkspaceCache(r"SDEConnection")
print arcpy.GetMessages() + "\n"

But this isnt working .. it crashes my Python window :(


texts: Getting signal when toggling group in PyQGIS?

texts:
I have a group which contains x number of layers. I can connect the group's node visibilityChanged signal to a simple print function so that when the group is toggled, a message is printed to the console.

The issue is that the message is printed x times (the number of layers in the group) + 1 (the group node itself).

Is there a signal which is only emitted once when the group node is toggled?

I tried searching in the QgsLayerTreeGroup and QgsLayerTreeNode classes but I seemed to have missed it.

A simple example:

def print_me():
    print 'something'

group = QgsProject.instance().layerTreeRoot().findGroup('myGroup')
group.visibilityChanged.connect(print_me)


texts: Algorithm split image doesn't show up anymore in QGIS 2.18.14

texts:
I developed a plugin for qgis in python language with a code line to execute split image algorithm from orfeo toolbox. When I update from QGIS 2.18.9 to 2.18.14 the split image algorithm doesn't show up in the processing toolbox and when I run the plugin I created, QGIS crashes because it can not find the algorithm.

I reinstalled QGIS and also reinstalled windows but I cannot find that algorithm in processing toolbox.

Is there an alternative to separate multispectral image in bands?


texts:
I am attempting to impute temperature based on elevation in ArcMap 10.2, using a regression equation. So far, I have been unsuccessful, as my output raster from my regression equation has the high too cool by 20 degrees, and the low too warm by 20 degrees.

I've attached my first output of the regression over the 30m DEM I've used

I've also attached the scatter plot of temperature points vs. altitude, which I then make the regression equation from. The equation I got from using a second order polynomial regression is as follows: -.00001*"DEM"^2+.0918*"DEM"

My approach is that there will be a relationship between decreasing temperatures and an increasing altitude, and that I can make a regression equation off of this.


texts: How to show average walking distance to closest point in census tract in ArcMap?

texts:
Given a map layer of point locations, and a second map layer of census tract polygons for the same area, how would you create a map layer that shows the average walking distance to the closest point for each census tract?

Euclidean distance?


texts: How do I get an accurate time attribute when importing a GPX layer as a vector

texts:
When I import a GPX file in QGIS (as a vector as track points) and I check the time attribute it appears to be 4 hours later than it should be. Interestingly when I take this same GPX file and use GPSBabel to convert it to a CSV file I find that the time listed in the CSV is the correct time.

So am I doing something wrong or is this a bug in QGIS.

Thanks for the help.

David


texts: ASC to JSON conversion. Where are the attributes?

texts:
I am converting MetOffice Nimrod .dat files to GeoJSON via the asc format.

The MetOffice dat format data is a proprietary type and I am using a script found here that converts from dat to asc.

Using QGIS I can see a reasonable looking output from this asc file.

From the dat file I use gdal to convert to json. I also do some geospatial clipping and convert to the epsg:4326 projection. I am confident from testing other files that these sections are correct.

I iterate over the asc files and use gdal, as is show below.

for file in glob.glob("*.asc"):
        print("Converting file " + file + " to json")
        src_ds = gdal.Open(file)
        if src_ds.RasterCount !=1:
            print "Raster issue - more than one band detected"
            sys.exit(1)
        if src_ds is None:
            print 'Unable to open %s' % src_filename
            sys.exit(1)
        try:
            srcband = src_ds.GetRasterBand(1)
        except RuntimeError, e:
            print 'Band ( %i ) not found' % band_num
            print e
            sys.exit(1)
        dst_layername = file[:-4]
        drv = ogr.GetDriverByName("GeoJSON")
        dst_ds = drv.CreateDataSource( dst_layername + ".json" )
        dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )
        gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

This outputs a reasonable looking GeoJSON file, as shown below. I am happy that this is geospatially correct.

However, my issue lies in the attributes for each of these polygons. The MetOffice NIMROD data contains precipitation measurements.

If I open the json file in a text editor, I see the following:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.1872638449, 51.653290137600003 ], [ -0.1728164805, 51.653066222 ], [ -0.1735389107, 51.635092885699997 ], [ -0.1879805698, 51.635316657899999 ], [ -0.1872638449, 51.653290137600003 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.1735389107, 51.635092885699997 ], [ -0.1590974392, 51.634867338200003 ], [ -0.1594612233, 51.625880718799998 ], [ -0.1738998445, 51.626106194099997 ], [ -0.1735389107, 51.635092885699997 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.1594612233, 51.625880718799998 ], [ -0.1450227909, 51.625653468899998 ], [ -0.1453892349, 51.616666906500001 ], [ -0.1598248185, 51.616894083799998 ], [ -0.1594612233, 51.625880718799998 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 0.1687182185, 51.539341131199997 ], [ 0.1831263429, 51.53907385 ], [ 0.1826964281, 51.5300889699 ], [ 0.168291136, 51.530356165599997 ], [ 0.1687182185, 51.539341131199997 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.4960296, 51.2999988 ], [ -0.4960296, 51.6806202 ], [ 0.275826, 51.6806202 ], [ 0.275826, 51.2999988 ], [ -0.4960296, 51.2999988 ] ] ] } }
]
}

There is no property or attribute metric here which gives a number to the precipitation for this polygon.

Where am I going wrong? Is there a further gdal argument for bringing the file attributes across from the asc file too?


texts: Determining length of roads that fall under low elevation using ArcGIS Desktop?

texts:
I am working on a project on sea level rise impacts of roads. Basically I need to determine length of roads that fall under low elevation (from 1m to 5m elevation).

I have a DEM file and road shapefile. I first converted DEM to integer and then raster to polygon. Then I tried to intersect my polygon with polyline but got error as output.

If I am going wrong somewhere can you suggest an alternative way of doing such analysis?


texts: Geoserver outputs latin1, even though utf-8 is set as encoding

texts:
I've been told, by the administrator, that this geoserver is outputting utf-8:
http://wfs-kbhkort.kk.dk/web/, however I can only get good results when parsing WFS from latin1 to utf-8. Can anyone confirm that the geoserver is correctly outputting utf-8 or confirm that it is latin1. If it is latin1, then how can I tell the administrator to change it to utf-8?

I know this perhaps is a rather odd question, but I'm very new to geoserver (and GIS in general)...


texts:
I'm trying to add a Layer Control (for streets and towns) in Leaflet but can't get it to work - not sure where I've gone wrong. My code is made up of two functions. One gets data using AJAX (JSON), and the other plots the data. I wont post all the code, just the parts relating to the layers:

function getJSON () {

  var map = L.map('map', {
  center: [57.73, -98.99],
  zoom: 10,
  layers: [streets, towns]
  });

  var baseMaps = {
  "Streets": streets
  };

  var overlayMaps = {
  "Towns": towns
  };

  var streets = L.tileLayer('https://api.mapbox.com/styles/v1/mapbox/streets-v10/tiles/256/{z}/{x}/{y}?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpejY4NXVycTA2emYycXBndHRqcmZ3N3gifQ.rJcFIG214AriISLbB6B5aw', {id: 'mapbox.streets',
      attribution: 'Map data © OpenStreetMap contributors, CC-BY-SA, Imagery © CloudMade'})

function plotData {

  var townMarkers = L.markerClusterGroup();

  var towns = L.layerGroup([marker]);

  for (var i = 0; i < townData.length; i++) {

  var townLocation = new L.LatLng(townData[i].lat, townData[i].lon);

  var marker = new L.Marker(townLocation);

  marker.bindPopup("<div class=townInfo><h1>" +
  townData[i].name)

  townMarkers.addLayer(marker);
  }
  map.addLayer(townMarkers);
  L.control.layers(baseMaps, overlayMaps).addTo(map);
  }
  }

I think I have one or two of the variables in the wrong location, but have tried them inside the for-loop and inside the plotData function.

Can anyone see where I have gone wrong?


texts: How to create an array of points in ArcMap?

texts:
I have a master point in following projection system:

WGS_1984_UTM_Zone_45N
WKID: 32645 Authority: EPSG

Projection: Transverse_Mercator
False_Easting: 500000.0
False_Northing: 0.0
Central_Meridian: 87.0
Scale_Factor: 0.9996
Latitude_Of_Origin: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
Spheroid: WGS_1984
Semimajor Axis: 6378137.0
Semiminor Axis: 6356752.314245179
Inverse Flattening: 298.257223563

Now i want to create 15 points toward north,15 points toward east,15 points toward west,15 points toward  south at 500m distance interval centering master point.These point has attribute to be populated in attribute table

Now how to create these points

Thanks


texts: Where is the Speed Unit option in Update Traffic Data tool of ArcGIS Network Analyst?

texts: I can't find the Speed Unit option in the Update Traffic Data tool for ArcGIS 10.2.0 (Enterprise). In the help Esri describes that you can specify a Speed Unit but I don't see the parameter.


texts: Add folder name as attribute to all shapefiles in that folder

texts:
I have many folders with shapefiles in each of them. I want to use python to go through each folder, add a new field to all the  shapefiles containing the folder name.

Any ideas?


texts: What should come first-Radiometric calibration or co-registration?

texts:
For a pair of images,in the pre-processing of satellite images, which step should come first-
Radiometric calibration or co-registration?


texts: Open source tool to convert a DEM to a 3D model

texts:
I would like to convert a DEM (Digitial Elevation Model) file (e.g. from SRTM) into a 3D model that I can then edit with regular 3D modeling tools (e.g. Meshlab/Blender/etc.). However I'm using Ubuntu Linux.

So, is there any open source programmes / tools that can do this? And if so, can you give me a simple guide for how to convert a DEM to a 3D model. (command line programmes perferred)


texts: How can I use the authagraph projection in QGIS?

texts:
I have only been able to find this resource which reproduces the authagraph projection and it is compiled in C.

Is there any way I can use this resource to make an authagraph projected map in QGIS?

Is there another way the authagraph can be made in QGIS?


texts: ArcGIS Server - host ArcGIS Services through fully-qualified domain name and servername

texts:
I’m in need of some direction…

We currently host an ArcGIS Server 9.3.1 SP2 .NET on our Intranet. A situation has come up that requires me to have the ArcGIS server respond to both the servername and through a fully qualified domain URL to access GIS services. How do I configure the server to respond to the fully-qualified domain name and servername for ArcGIS requests?

http://servername/ArcGIS/rest/Services

http://sub.domain.foo/ArcGIS/rest/services

Is there a whitepaper, directions, something that explains how I can fulfill this scenario?

mwalker - Thank you for the response.

Reverse proxy would be ideal and actually after reading through this setup I will be proposing we setup our environment as such going forward.

In the mean time I need to have this particular Silverlight site up and running by Monday morning within our Intranet for remote office access (Secure VPN), using a single domain name where the GIS services can be referenced through the domain name itself.

Now I know I could bind the domain name to the Default Web Site within IIS, setup the Silverlight project, and GIS services and the site should respond as expected. This is the quick and dirty answer to my dilemma for a Monday morning availability.

But for conversation and learning more about the ArcGIS server flexibility - are there any other options to consider?


texts: How to buffer increasing/decreasing population zones?

texts: I will find out the difference in population between year 2000 and 2013. I have a feature class file with the population in year 2000 and 2013 that are connected to the municipal number and that is a point theme. What I want is to have buffer zones around each point and the size of the buffer zone to be determined by the number of residents that each point has. If the population has increased the buffer zone will get a red color and where the population has decreased the buffer zone will be blue. Is there anyone who can help me with this? Whether as a model in ModelBuilder, Python or a Toolbox.


texts: Converting Number to Percentage in ArcGIS Field Calculator?

texts:
Using ArcGIS 10.2 and Personal geodatabase,I have a field with values of Numbers like  0.461538 and Type of "Double".

How can I Export them to new Field in Percentage Format to be shown like 46% or at least round the number to looks like 46 then I can add "%" on label to them?

(Honestly! I do not think this is a good idea for further Analysis!).


texts: How to show/hide labels with mapserver and openlayers

texts:
I'd like to have a check/uncheck option like this print checkbox to show/hide labels in my WMS layers served by mapserver 5.6.7 in openlayers 2.11. I have a toolbar where to put this using probably geoext1 or ext3.4, so the questions:

Is this possible?
If so, how can it be done?


texts:
I have created a map with multiple points and path layers in QGIS 3.8 (Zanzibar). When I try to publish the map with the QGIS2Web plugin, there are two kind of results:

When I use the "OpenLayers" option, it works quickly and it shows the raster map and the paths perfectly; however, only four of the many (>100) points are shown.
When I use the "Leaflet" option, it does show all the points (although extremely slowly), but no raster map, and the functionality is really poor (very slow and barely visible).

I prefer to use the OpenLayers option, as it worked perfectly, except for all the missing points and labels.

Could anybody tell me how I can fix this?

The underlying dataset seems to be fine, as all the points are shown in the Leaflet option, and also when I load them in Geojson.


texts: Converting local solar time to UTC for MODIS

texts:
I am having a hard time to convert local solar time to UTC for MODIS data.

According to the Users' Guide of the product (MOD11A1 Product - Wan Z., 2013), local solar time given for each pixel is

  "the UTC time plus grid’s longitude in degrees / 15 degrees (in hours,
  +24 if local solar time < 0 or -24 if local solar time >= 24)".

Should we take into account the day of year when calculating the UTC time although it is not mentioned in the instruction provided above?

I am using Google Earth Engine, code editor with Javascript.


texts: Error: already exists.(stg_e_filealreadyexists) when using IMapDocument.Open method?

texts:
I got the error: already exists.(stg_e_filealreadyexists) when using the method IMapDocument.Open() as below codes:

if (File.Exists (strDirectory + "\\temp.mxd"))
    File.Delete(strDirectory + "\\temp.mxd");
System.IO.FileStream oFileStream = new System.IO.FileStream(strDirectory + "\\temp.mxd", FileMode.Create);
oFileStream.Write(MyImage, 0, MyImage.Length);
IMapDocument mapDoc = new MapDocumentClass();
mapDoc.Open(strDirectory + "\\temp.mxd"); // Throw exception

Do you know why and how to fix this error ?


texts: Why do the rasterio Resampling methods ‘max’, ‘min’, ‘med’, ‘q1’, ‘q3’ give a segmentation fault?

texts:
I want to resample a georeferend raster .tif file with Python using rasterio. Here is the example I found example and the documentaion of the resampling methods.

When I am using the following methods ‘max’, ‘min’, ‘med’, ‘q1’, ‘q3’ which are only supported in GDAL >= 2.0.0 the program fails.

I am using a conda environment and rasterio.gdal_version gives me version 3.0.4.

When I run the following code (i.e. almost the example code):

raster = Path('dsm_EPSG_4839_cropped.tif')

import rasterio
from rasterio.enums import Resampling

upscale_factor = 1/(2**5)

with rasterio.open(raster) as dataset:

    # resample data to target shape
    data = dataset.read(
        out_shape=(
            dataset.count,
            int(dataset.height * upscale_factor),
            int(dataset.width * upscale_factor)
        ),
        resampling=Resampling.max
    )

    # scale image transform
    transform = dataset.transform * dataset.transform.scale(
        (dataset.width / data.shape[-1]),
        (dataset.height / data.shape[-2])
    )

I am getting the following error:

Segmentation fault: 11

How do I solve the issue and get the example to run without the segmentation fault?


texts:
I followed the instructions manually building a tile server and ran in some problems.

When entering http://localhost:8080/osm_tiles/0/0/0.png

renderd shows this:

renderd[263]: DEBUG: Got incoming connection, fd 22, number 9
renderd[263]: DEBUG: Got incoming request with protocol version 2
renderd[263]: DEBUG: Got command RenderPrio fd(22) xml(default), z(0), x(0), y(0), mime(image/png), options()
renderd[263]: DEBUG: Failed to read cmd on fd 22
renderd[263]: DEBUG: Connection 8, fd 22 closed, now 8 left

It seems renderd is aware of the request, but there do never appear some files in /var/lib/mod_tile. The error.log of the apache webserver shows:

debug: init_storage_backend: initialising file storage backend at: /var/lib/mod_tile
[Mon Jun 01 12:33:15.828298 2015] [tile:info] [pid 104:tid 140230287886080] [client 172.17.42.1:60074] tile_storage_hook: handler(tile_serve), uri(/osm_tiles/0/0/0.png)
[Mon Jun 01 12:33:15.828355 2015] [tile:info] [pid 104:tid 140230287886080] [client 172.17.42.1:60074] Requesting style(default) z(0) x(0) y(0) from renderer with priority 5`

and then responds me with a 404 Page not found.


texts: Getting results from processing.run GRASS r.resample in PyQGIS

texts:
I am trying to get output from r.resample in a stand-alone PYQGIS application. The code runs and does not show any errors, but it does not save an output file (i.e. 'c:/path/to/output.tif'). As far as I can tell, I have set GRASS up ok. I can run the same code through the QGIS console and get a result.

I found a possible answer here, saying that it might be because grass commands need an extent. However, I tried that (see code below) and the command still does not run.

I am using Windows 10, Python 3.7 and QGIS 3.6.

What is the reason it is not producing an output raster?

from sys import argv
import os, sys, time, csv, ogr,subprocess
from osgeo import ogr

# GRASS7 set-up ################################
gisdb = os.path.join(os.path.expanduser("~"), "grassdata")
if sys.platform.startswith('win'):
    # MS Windows
    grass7bin = r'C:\OSGeo4W64\bin\grass76.bat'
# query GRASS GIS itself for its GISBASE
startcmd = [grass7bin, '--config', 'path']
try:
    p = subprocess.Popen(startcmd, shell=False,
                         stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out, err = p.communicate()
except OSError as error:
    sys.exit("ERROR: Cannot find GRASS GIS start script"
             " {cmd}: {error}".format(cmd=startcmd[0], error=error))
if p.returncode != 0:
    sys.exit("ERROR: Issues running GRASS GIS start script"
             " {cmd}: {error}"
             .format(cmd=' '.join(startcmd), error=err))
out2 = out.decode("utf-8")

gisbase = out2.strip(os.linesep)
print("GIS base: ",gisbase)
os.environ['GISBASE'] = gisbase
# define GRASS-Python environment
grass_pydir = os.path.join(gisbase, "etc", "python")
sys.path.append(grass_pydir)
# GRASS7 set-up complete #########################

app = QgsApplication([], False)   # second argument to False disables the GUI
app.initQgis()
# start processing
import processing
from processing import Processing
Processing.initialize()

#### r.resample starts here ####################
raster_source = path/to/raster.tif
rlay = QgsRasterLayer(raster_source, "raster_name", "gdal")
if not rlay.isValid():
    print("Layer failed to load!") # layer is valid

ext = rlay.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
coords = "%f,%f,%f,%f" %(xmin, xmax, ymin, ymax)

# You can see what parameters are needed by the algorithm
# using: processing.algorithmHelp("grass7:r.resample")
 params = {
     'input' : rlay,
     'output' : 'c:/path/to/output.tif',
     'GRASS_REGION_PARAMETER':coords,
     'GRASS_REGION_CELLSIZE_PARAMETER':200,
     'GRASS_RASTER_FORMAT_OPT':'',
     'GRASS_RASTER_FORMAT_META':''
                    }

res = processing.run("grass7:r.resample", params)
## Process runs but no output is saved

Edit: some commands output the final file, but not all do. For example, gdal:proximity also does not provide an output.


texts: Counting traversed arcs on routes to write attributes using ArcGIS for Desktop?

texts:
I have number of persons located on a road network and I have calculated their shortest distance to a particular location.

From the routes output created from Closest Facility can I

break the routes down into the individual arcs traversed by each person?
assign each arc a count attribute indicating the number of individuals who use that arc to get to their destination?

I am using Network Analyst in ArcGIS 10.2 for Desktop.


texts: Assigning QgsCategorizedSymbolRenderer Spectral Ramp to MultiPolygon

texts:
I currently have a Python script that adds a MultiPolygon shapefile to my map within QGIS. Within QGIS, I can apply a categorized symbology render, with the color ramp 'Spectral' going from blue to red, with 55% opacity, for a shapefile called 'Shapefile'. The values within the file under the field 'Intensity' are applied within the categorized symbology (see screenshots). However, I would like to apply this within the Python console. Is there any way I can do this?

UPDATE: Using the code suggested below, I have this:


texts: How to find most common classification within buffers?

texts:
I have a point shapefile that identifies the location of incinerators in England and an output area classification shapefile.

I have created a buffer for all the points in the shapefile because i want to see what area classifications fall within the buffer of all the incinerators i have mapped.

Next, I intersected the buffer i have created with the output area classification layer.

What I want to know is how to i analyse the data within the buffers of the incinerators?

I want to to know what is the most common classification throughout all the buffers....


texts: Convert latitude longitude pair to pixels in geotiff

texts:
I am not a GIS guy.

I have a GeoTiff and I need to read elevation information out of it using GDAL (in some Java Code). So now I need to convert a given latitude/longitude pair into a pixel inside the GeoTiff.

Here's some information gdalinfo gave me about the file:

Driver: GTiff/GeoTIFF
Size is 58808, 30323
Coordinate System is:
PROJCS["MGI_Austria_Lambert",
    GEOGCS["GCS_MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel_1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",46],
    PARAMETER["standard_parallel_2",49],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",400000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (106549.267203768889885,576922.512073625810444)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  106549.267,  576922.512) (  9d19' 7.63"E, 49d 1'24.32"N)
Lower Left  (  106549.267,  273692.512) (  9d31'19.50"E, 46d17'55.02"N)
Upper Right (  694629.267,  576922.512) ( 17d21'50.31"E, 49d 1'22.34"N)
Lower Right (  694629.267,  273692.512) ( 17d 9'35.51"E, 46d17'53.15"N)
Center      (  400589.267,  425307.512) ( 13d20'28.29"E, 47d43'39.80"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
NoData Value=-3.40282306073709653e+38

Can somebody tell me the formula on how to convert my coordinates into pixels?


texts: Raster smoothing - assigning each raster the average value of its neighbors

texts:
I have a raster world map of population densities at a resolution of 30 arc seconds. For each raster, I want to find out the average density of neighbors within a certain distance (e.g. 0.1 degrees). The aim is to keep and save in an extra file only those raster values with a density greater than a certain threshold (e.g. 1000 people per sqkm). All contiguous rasters above this threshold will then be clustered to form a core city.
Next, I want to find the city fringe to those core cities. This shall be composed by rasters (attached to the core cities) with population densities between 500 and 1000. Core plus fringe will compose a complete city. Those cities then should be transformed to a polygon shape file where I can assign each city an own ID. Furthermore, having a polygon shape file allows me to easily assign values from other raster maps to the city using zonal statistics.

My question is how to do all of this. I tried around a little with Warp but I think that the result is different from what I want because Warp only is able to average over the neighbors to reproject the small to bigger raster, but which rasters are defined as neighbors depends on the layer extent chosen.

Below a sample image of my raster map.

Is there a way without using GRASS? Because I cannot get GRASS running..


texts: ArcGIS Server - writing error messages to jobs

texts: I'm trying to aid some of my other teams by improving the error handling when geoprocessing services are being called. In my code I'm handling errors however after publishing that code to the server 10.8. I cannot see any errors being written out under jobs\jobid URL of the REST API. Is there some specific standard I need to follow in order for messages to be written to ArcGIS Server REST API when looking at the status of a job?


texts: Removing a plugin from custom menu in QGIS3

texts:
I am developing a plugin in QGIS3 and I want it to be located on it's own menu. I am able to create the plugin and load it just fine.  But when I try to unload the plugin, nothing happens. Then if I turn the plugin on again a second entry appears in my menu.
How do I unload a plugin on a custom menu in QGIS3?
def initGui(self):
    # Check if the menu exists and get it
    self.menu = self.iface.mainWindow().findChild( QMenu, "&Station Tools" )

    # If the menu does not exist, create it!
    if not self.menu:
        self.menu = QMenu("&Station Tools", self.iface.mainWindow().menuBar())
        actions = self.iface.mainWindow().menuBar().actions()
        self.menu.setObjectName( '&Station Tools' )
        lastAction = actions[-1]
        self.iface.mainWindow().menuBar().insertMenu(lastAction, self.menu)

    self.action = QAction(QIcon(":/plugins/delete_feature/icon.png"),"&Delete Selected Feature", self.iface.mainWindow())
    self.action.triggered.connect(self.run)
    self.menu.addAction(self.action)

    # will be set False in run()
    self.first_start = True

def unload(self):
    """Removes the plugin menu item and icon from QGIS GUI."""
    for action in self.actions:
        self.iface.removePluginMenu('&Delete Selected Feature',action)
        self.iface.removeToolBarIcon(action)

I believe my initGui works but it is my unload where I am not doing something right because the tool will not unload from the QIGS UI.


texts: How to use OSM xml extracts? Can I make an offline app in Android using OSM?

texts:
I'm new to OpenStreetMaps. I have followed their wiki/main page and tried to understand it.

Can anyone briefly explain how exactly OSM xml data works and how it is rendered? ...Any links /pointers would be useful.

I went to this link: cloudmade osm xml extracts.  Here I'm seeing various file extensions of the OpenStreetMap files.  I want to know, is there a way to use OSM xml extracts to download /display a map?

Also, how can I use OpenStreetMap data to make an Android app completely offline?
E.g., tripadvisor recently launched such apps.

Any help, suggestions are welcome.  Thanks.


texts:
First of all, I really could not think of a good title for this question, so feel free to modify it.

I keep encountering an error using ArcMap 10.5.1. I'll have a table open, to which I add a field.

Here is a table to which I added the fields "x" and "y" for the x and y coordinates:

Then when I attempt to "Calculate Geometry" for either of these new fields, I get the following sequence of error messages:

I also tried to use the "Add Surface Information" tool in 3D Analyst, and it gives me the following error messages:

I also tried adding a field to a polygon shapefile and then adding surface information. It gave me the following error message, which may actually point to the source of the problem. It appears to be a mismatch between the number of shapes vs the number of records. I have no idea why, since the attribute table appears to look correct:

One curious thing is that the shapefiles appear to be corrupted after this occurs. Does anyone have any ideas?


texts: Is it possible to cut/filter all layers to an area of interest in one step?

texts:
The problem of working with a lot of data (about 100 layers). Very slow computer, rendering long.

Vector -> Geoprocessing Tools -> Clip a very good tool, but replace the 100 layers is very high and can go wrong. In addition to work with the same data in different regions, and each time the clip function is complicated work.

Is there a faster way, and saves computer resources?

On ArcView GIS 3.2 is called the "Area of interest" or Autocad Map 3D "query to filter date".


texts: Use try, catch when querytask times out using JS API

texts:
When I do a querytask and the date range criteria used is too vast, I get a 500 internal server error after 1.7 minutes every time (unless the query takes less than 1.7 minutes, then it runs fine). Is there any way I can use try, catch to let the user know the search has timed out.

I tried the following try, catch code but it does not do anything. It still says "Searching..." in my 'wait' div.

    function execute () {
            query.geometry = addressPoint;
            startDateQuery = document.getElementById("startDate").value
            endDateQuery = document.getElementById("endDate").value

            //esriConfig.defaults.io.timeout = 60,000;
            if (endDateQuery !== null && endDateQuery !== undefined && endDateQuery !== ""){
                query.where = "(DATE_VAL >= '" + startDateQuery + "') AND (DATE_VAL <= '" + endDateQuery + "')";
            }
            else {
                query.where = "DATE_VAL = '" + startDateQuery + "'";
            }
            try{
                queryTask.execute(query, showResults);
                document.getElementById('wait').innerHTML = "Searching...";
            }
            catch(err) {
                document.getElementById('wait').innerHTML = "Search timed out. Use a smaller search range.";
            }
        }


texts: Determinate Parcel Neighbors using PYTHON

texts:
As you can see in the First picture, I have parcel 220, and i want to determinate the position of each of his neighbor ! So the idea was to calculate the bearing between the centroinde of 220, and the other centroid of neighbors (when the bearing is between [50,150] = North; [150, 250]= East; [250,350]= West ...). So to determinate the neighbors i use the function arcpy.PolygonNeighbors_analysis (in_features, out_table, {in_fields})(Picture 2), and at the picture 3, we have a parcels with centroide (X,Y).
So each time i should take

(X,Y) of P220 as start Point and (X,Y) of P228 as end point than calculate bearing
(X,Y) of P220 as start Point and (X,Y) of P227 as end point, calculate bearing
(X,Y) of P220 as start Point and (X,Y) of P226 as end point, calculate bearing
than for all the neighbors on the Picture 2

For doing this I used three loop for, and I think it's turning right, but it does it just for the first row, and doesn't move to the newt row!

The code is this:

>>> cur = arcpy.UpdateCursor("Join_Output_12")
>>> cur2 = arcpy.UpdateCursor("Join_Output_12")
>>> cur_t = arcpy.UpdateCursor("tb20_12_2014_")
>>> del row2, Refname__table_neighbor, start_x1, start_y1, end_x1, end_y1, copy_field
>>> for row in cur :
...     Refname_feature = row.getValue("RefName")
...     for row_t in cur_t :
...         Refname__table_target = row_t.getValue("src_RefName")
...         if(Refname__table_target==Refname_feature):
...             start_x1 = row.getValue("POINT_X")
...             start_y1 = row.getValue("POINT_Y")
...             print Refname_feature
...             print start_x1
...             print start_y1
...             Refname__table_neighbor = row_t.getValue("nbr_RefName")
...             print Refname__table_neighbor
...             for row2 in cur2 :
...                 copy_field = row2.getValue("copy_field")
...                 if (copy_field==Refname__table_neighbor):
...                     end_x1 = row2.getValue("POINT_X")
...                     end_y1 = row2.getValue("POINT_Y")
...                     print copy_field
...                     print end_x1
...                     print end_y1
...     row = cur.next()

The result is:

1
221408.8543
200288.7235
1
1
221408.8543
200288.7235

So I think that the problem comes from for loop!

Now after getting this result, that you can see on the picture 4 (Position off each parcel)

The result i get was good! But it doesn't respond at 100% to my need, to be more clear let's see the example of parcel 221, it has parcel neighbor 220 at North and also West, but at the table we can see just North, for that i want to canncatenate this result to the result of the function CalculateAdjacentFields_cartography (in_features, in_field)
you will tell me why just not to use only this function ? but it's also uncomplete ! look at picture 5 to see the result

Let's take the exmaple P220, we can see that it sets just 224, 226 on North, even if there is More !
So the idea is to concatenade the two results to get at last one more complete and reliable result !
But how to do it ? i can't find one logic to combinate them !!

EDIT :
This is an example of complicated relationship, and on the table after what i get as result

the lin in color, presents Neighbors picked up by the program, and when i looked what a kind of parcels mess, i found is just who is not in touch with the taget parcel (as 164 P2 and 150 in North..)!!


texts: Is the pixel aspect ratio of a raster dataset always 1:1?

texts: When we are talking about rasters, is the spatial resolution of the pixel defined in a raster dataset always fixed to a ratio of 1:1, eg, 1m x 1m ?


texts: GeoWebCache returns files as geoserver-dispatch.application instead of kml in Chrome

texts:
This url: http://maps.gcc.tas.gov.au/geoserver/gwc/service/kml/GCC_cuo:GCC_Image.png.kml should download a kml file, but it doesn't in Chrome. It downloads a geoserver-dispatch.application file instead.

It does download fine in IE.

Is there something I can do to fix this?


texts: Show/hide individual selected feature qgis

texts:
Is it possible to select individual features of a layer on a map and hide them? (I am working with
3.8.3-Zanzibar version of QGIS).

I need to print different maps of a close by area and I want to show only the symbols that are relevant in each case, but not others close by.

In the example below, I have selected the points (in yellow) that I would like to momentarily hide.

I know I could add a column in the associated table with a numeric field and fill in the ones I need to show with 1, the rest with 0 and then apply a filter. But I was looking for something more direct, like the show/hide label button.


texts: Verifying formula that will convert longitude 0-360 to -180 to 180?

texts:
I have a netCDF file that contains climate data in a Gaussian grid. The longitudes are given from 0 to 360. I need to convert this to -180 to 180.

Is this the correct formula -

If longitude is greater than 180 subtract 360

I am looking to do this in Fortran.


texts: How to get an indication of how points are spatially distributed?

texts:
Let's say I have a set of points and I want to get an indication of their spatial distribution. That is, if all points are at one location a possible return value would be 0, whereas if all points are distributed evenly the return value would be 1.

The use case: I have two data set which cover the same area. I'm developing a tool to compare both sets with each other. I also want to include a rough indication of the point distribution to detect differences in the sets.

My program will be written in python using open source tools (so no arcgis).


texts: Operations Dashboard - How to insert a line break into summary widget trailing text?

texts:
In a ArcGIS Online Ops Dashboard Operation View, I have a summary widget configured to show a count of features like so:

6,087 new features added

...where 6,087 is the count of my features and " new features added" is my trailing text that I add into the configuration.

What I would like it to look like is:

6,087 new features added

In other words, I want a line break in it. I've tried adding these into my trailing string in the widget settings:

\n
<br>
&lt;br&gt;
%26lt;br%26gt;

...and they all just come out raw. Also, I've tried adding in the breaks in both the Ops Dash GUI app and I edited the Operation View JSON directly in AGO Assistant. Anyone got a trick to get a line break in there? Not even sure if it's possible.

A colleague jogged my brain and got me to look at the HTML output for that element, which looks like so:

<text dominant-baseline="middle" class="label" fill="#f79646" text-anchor="middle" x="189.5" style="font-size: 20px;">6,087 total %26lt;br%26gt;dead meters</text>

So, it appears that my text all gets wrapped up in the <text> tag, which would explain why none of my HTML takes. I'm beginning to think this would be more work than it's worth to accomplish.


texts: Migrating ArcGIS Server 10.02 to a new server

texts: Our server that runs ArcGIS Server had major issues this morning, sending us into a panic to get it installed on a new server (as was in our plans in the coming weeks). Now that we are ready to move all our services/applications over I've come to realize I have no idea where to begin. Is there any method to copy the data over to the new server and still be able to use ArcGIS Server without having to re-create all our applications? I copied the services and applications, but, not to my surprise, it didn't work. I tried IIS Web Deploy, but run into issues. Any help would be greatly appreciated!


texts:
I am trying to convert shapefile that was created by FME in EPSG:900913 to EPSG:3857 (I know its the same definition) with OGR in a script.

datasource = ogr.Open(os.path.join(dirname,shapefileName))
layer = datasource.GetLayer(0)
srcSpatialRef = layer.GetSpatialRef()
dstSpatialRef = osr.SpatialReference()
dstSpatialRef.ImportFromEPSG(3857)
coordTransform = osr.CoordinateTransformation(srcSpatialRef, dstSpatialRef)

for i in range(layer.GetFeatureCount()):
    srcFeature = layer.GetFeature(i)
    srcGeometry = srcFeature.GetGeometryRef()
    srcGeometry.Transform(coordTransform)
    geometry = GEOSGeometry(srcGeometry.ExportToWkt())

I got this error:

  GDAL_ERROR 6: No translation for Mercator_Auxiliary_Sphere to PROJ.4
  format is known.

This is weird because the scripts works if the shapefile is not created by FME but by QGIS both in 900913.

this is the .prj of the not working shapefile

  PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG",3857]]

Does it means that OGR cannot transform file from FME?


texts: Forcing polyline to offset from polygon border

texts:
I am currently working on a site with polygons and polylines and am looking for a way to ensure that the polyline does not overlay the site boundary.

I am creating a habitat map and do not want the hedges/ditches which are present on the edge of the site from obscuring the site boundary or vice versa. This only affects some of the polylines, as obviously those in the centre of the site are not affected. Moving the lines manually seems clunky, time-consuming and also wouldn't work at certain zoom levels. Any ideas?


texts: QGIS time manager: export video on Windows

texts:
Recently, the time-manager-pugin seems to support video-export on Windows: https://github.com/anitagraser/TimeManager
However, I can't find any instructions how to install this feature. I updated to time manger plugin, however, it still mentions "only Linux/OSX" - so how to install or activate this feature?


texts: Looking for web service to get State, County, and Place FIPS codes from an address

texts:
I have been searching for a few days to find a web service that will provide the State, County, and Place FIPS codes from an address input. Most of the services I've found only provide a State, County, and Census block code.

For example, for the address 4139 S 143rd Cir, Omaha, NE 68137 would return:

State: 31
County: 05
 Place: 537000


texts: Map scale vs. map extent vs. map canvas size in Print Map Composer

texts:
I am using QGIS 2.0 and when I try to create my map in the print composer, I cannot make the scale and the extent the way that I need it.

I need my map to be at the scale of 1:10 000. However, the map area is quite large and I need a specific x min and y max in the extent (the x max and y min do not matter).  I also need a full A4 size sheet of paper.

When I enter the scale in, the map extent is not where I need it.  If I change the extent, it either changes the scale of the map or changes the size of the map canvas in the print composer.  When I change either one of those back, the scale changes and so on and so forth.

Is there some way to lock the scale and then just move the map image inside the map canvas?


texts: OSM Map projection doesn't seem to change

texts:
I have the following OpenLayers map:

Proj4js.defs["EPSG:31258"] = "+proj=tmerc +lat_0=0 +lon_0=13.33333333333333 +k=1 +x_0=450000 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs";
var googleProj = new OpenLayers.Projection("EPSG:900913");
var salzburgProj = new OpenLayers.Projection("EPSG:31258");
var bounds = new OpenLayers.Bounds(426249,294999,426251,295001);
var maxZoom = 19;
var minZoom = 15;
var res = [4.77731426696777,2.38865713348389,1.194328566741945,0.5971642833709725];

mapCOD = new OpenLayers.Map({
    div: "mapCODs",
    allOverlays: true,
    maxExtent: bounds,
    projection: salzburgProj,
    allowSelection: true,
    controls: [new OpenLayers.Control.PanZoomBar(), new OpenLayers.Control.Navigation()]
});

var layer = new OpenLayers.Layer.OSM(
        "layer1",
        "tiles/input/${z}/${x}/${y}.png",
        {zoomLevels: maxZoom, zoomOffset: minZoom, resolutions: res}
    );
mapCOD.addLayer(layer);

I have produced tiles with mapnik using EPSG:31258, so I would like to change the default projection of the map from 900913 to 31258.
The problem is that "projection" parameter is not taken into account and if I perform a console.log(mapCOD.getProjection()), I obtain "EPSG:900913".

Another proof that map projection remains 900913 is that I have to perform a coversion to get the correct map center:

var center = bounds.getCenterLonLat();
mapCOD.setCenter(center.transform(salzburgProj,googleProj), minZoom);

What am I missing?


texts: Problem with creating a "blank" basemap in an application built with arcgis javascript api

texts:
I'm building an application where I'd like the user to be able to select a blank basemap from a gallery of basemap options.

I've found that if I initialize the app with a "white" "blank" or "none," I get the blank basemap I desire.  Unfortunately, when I set this property it seems to "lock" the basemap into the "blank" property. When the user selects a new basemap, such as "topo" from the gallery I've created, the basemap stays blank.

On the other hand, if I initialize the app with a basemap such as "topo," the user can switch between other basemaps, but if they try to select "blank" the basemap doesn't change.

Even though it isn't listed as an option in the API, is it possible to switch between a "blank" basemap and the others?


texts: How to convert Esri Compact Cache to other format? (mbtiles, et.c.)

texts:
I need to convert Esri Compact Cache-files to a "better" (a more open & common) format. Ideally I would like to convert to mbtiles or just 'tiles on disk'.

What software can perform this conversion, and how is it done?

The Compact Cache-files I'm working with can be downloaded from this site (see "Offline products"):
https://www.basemap.at/index_en.html#sec-zugang


texts: Raster layer - changing values in raster cells

texts:
I have a lot of cells in a raster layer (GeoTIFF) with a -9999 value in it. I would like to update that values to 0.

I mean, in vectorial layers would be something like: CASE WHEN band1=-9999 THEN 0 ELSE band1 END

Is this possible to do for raster layers in QGIS or any other GIS software?


texts:
I compiled GRASS 7 on Ubuntu 12.10, with liblas in order to run r.in.lidar and v.in.lidar. I now want to import .las data; the option in file -> import raster data -> LAS LIDAR points import is given and enabled, but when I select it nothing happens.

It seems to be a problem with libLAS. The console error message says:

  Unable to fetch interface description for command 'r.in.lidar'. Details: r.in.lidar: error while loading shared libraries: liblas_c.so.2.2.0: cannot open shared object file: No such file or directory.

Addition: running lasinfo from terminal works.

Any idea where the bug could be?


texts: Is there a limit to the number of raster layers in Raster Calculator?

texts:
I have 16-day Enhanced Vegetation Index (EVI) images for an area across roughly 10 years. 244 images in all. I would like to average these over the time period to get an average EVI value.

Is there a limit to the number of rasters ArcGIS Desktop can perform analysis on in Raster Calculator?


texts:
I am looking for solution to style the symbol in MapServer; with polygon outlined with simple line along with filled triangle inward like below

; I tried creating symbol with following way
SYMBOL
  NAME "Inward line"
  TYPE vector
  POINTS
   0   0
    5  10
    10  0

  END # POINTS
  ANCHORPOINT 0.5 0
END # SYMBOL

with line symbols adding together,. but no luck,.


texts: How to create a hole in a polygon in QGIS?

texts:
I'm trying to use the tool "delete ring" in order to create a hole in a polygon. The geometry of the hole is another feature of the same shapefile

I cannot manage to delete the part of the polygone. Does someone have a tutorial or some explanation?


texts: Is it possible to add house numbers using the Mapbox Street Vector Tiles source in Tilemill2?

texts: Title says it all. Does anyony know if the house number is included in Mapbox’ vector tiles, and if so how to access it for styling? Thanks for any advice!


texts:
I'm drawing many blocks like these

Then I need to surround them by points at fixed distances from each other like so

But I have a lot of blocks, and many points to calculate the distance, and then draw. I normally use the measure tool around the block, and then add the points on top

This process in extremely tedious, so I had a thought of drawing a line and by a certain way break that line into points that are separated from each other by certain distances and transferring those points into a different layer.

Is there anyway this way can work or does anyone have a better idea?


texts: Weighted zonal statistics in Python and/or R

texts:
I have raster files with precipitation data (4x4 km), monthly, for about 100 years. I have a shapefile for a watershed (which includes several polygons as sub-basins).
I need the average precipitation for the watershed for each month. I don't want to do this by hand in ArcMap because it will probably need to be repeated plenty of times (for other watersheds, and when more data becomes available), so I am coding it.

I have written code in R now, using packages raster and rgdal. It works great in the sense that it has a functionality to calculate the weighted average per sub-basin; it checks which fraction of a raster cell falls within the polygon and uses this to calculate a weighted average for the polygon. Because some of the sub-basins are small in comparison to the raster, this is a great feature.
And in addition, because I can access the shapefile information such as the area of the sub-basins, I can calculate a weighted average for the whole watershed.

However - this is slow. It also says it in the description of the function extract(). Calculating the weighted value takes 4-5 seconds, so when I ran this script for 30 years worth of monthly data, it took 37 minutes (there are a few more calculations besides the weighted average).

I have a feeling it might be faster to use Python, but I am not sure how to approach this. I did get started with this before the R code (see below) but shifted my attention to R. It seemed a bit daunting since it's been a while that I coded in Python, and I've never used Python for geospatial analyses.

Also, an initial test seems to show that in Python, there is no weighted averaging happening, and that small polygons that do not contain the center of a raster cell, get assigned 0 or None.

So my questions are:

If it is possible in Python, will weighted averaging that way be faster?
How to access other data associated with the shapefile (such as area)?

Below is the very very basic code I have so far.

    from rasterstats import zonal_stats
    import gdal

    basin="Steinhatchee_Project.shp"
    rain="PRISM_ppt_stable_4kmM3_201504_bil.bil"

    # NOTE From Python documentation: "There is no right or wrong way to rasterize a vector.
    # The default strategy is to include all pixels along the line render path (for lines),
    # or cells where the center point is within the polygon (for polygons)
    stats=zonal_stats(basin,rain)

    # get average for every polygon (there are 19)
    data=([f['mean'] for f in stats])

    # take out the "None" values...
    data2=[x for x in data if x is not None]

    # ... so I can calculate the mean for the whole watershed
    mean(data2)

data looks as follows. Note the None...

    Out[76]:
    [115.64999389648438,
     83.07333374023438,
     92.33000183105469,
     73.35416666666667,
     97.36499786376953,
     92.03500366210938,
     80.30999755859375,
     102.8699951171875,
     94.3499984741211,
     110.99250030517578,
     80.18000030517578,
     67.6749979654948,
     107.825,
     82.63999938964844,
     88.375,
     59.769999186197914,
     63.48249816894531,
     None,
     None]

My system specs: Windows 7, Python 2.7.9 (using the Enthought Canopy IDE) 32-bit, GDAL-1.11.3, rasterio-0.29.0, Shapely-1.5.13


texts:
I have a friend with a bit of a problem. He has 1000+ shapefiles mapping paddocks and being used for precision agriculture. He's recently discovered the base station positions used to create these maps were incorrect.

The base stations are being updated, which means each shapefile needs an x and y offset applied to all features. The offsets vary but are typically small (~1-2 meters). The correction needs to be applied with sub-centimetre accuracy.

I've been brought on board to automate the process. I've written a C# application that can accurately calculate the offsets to apply to each file. What I'm looking for is something to apply these offsets for me. Either a C# library or an external command line application would be ideal. I could implement something myself but I don't want to re-invent the wheel.

Does anyone have a suggestion for achieving this?

I've read that I can use ogr2ogr and specify a different false easting/northing for the source SRS, although I'm not sure if this will provide the accuracy I require. I also don't know what projection the files are using.

Any help would be greatly appreciated!


texts:
Im loading some layers from postgis using pyqgis in my standalone qgis app. For each layer i also load a style using:

lyr.loadNamedStyle("style.qml")

I've setup 2 versions of styles in qgis desktop, one with and one without scale range visibility settings for the layer. The problem is, when i use the style with scale visibility (no matter if i set a really large range from 1:100 000 000 to 1:1) my features are not shown, not ever. If i use no style at all, or style without scale visibility settings, the layers are visible normaly.

I should add, styles work in QGIS desktop as expected.

What bothers me slightly is that when i print the current scale property of the map canvas in my app, the value has an insanely large value, for example 10e14 to begin with, even though the scale should be somewhere around 4000000 (1:4000000). What is even wirder is, this scale value (canvas.scale()) changes significantly when i pan (not zoom) the map.

Any ideas what is going on?

Because i thought the mapCanvas srid is not setup correctly and the scale would never fit into my scale range, I tried setting map canvas srid explicitly, did that on layers, stil no luck. I am using QGIS 2.8.1

UPDATE: It seems in any case i have to set the CRS explicitly

canvas.mapRenderer().setDestinationCrs(QgsCoordinateReferenceSystem(32632))

but it doesn't change anything. I just thought that layers will implicitly enforce my local CRS but they don't.


texts: Clean up pixel noise from raster in ArcGIS

texts: I have this raster. It has some pixels near the main objects, this raster is a product of a classification. those little pixels have the same value of other big objects, so I cannot just reclassify. Is there any tool to do this clean up on this raster?.


texts: Saving active view into ArcGIS Pro project using ArcPy?

texts:
I have a desirable requirement to be able to set the active view (Layout, Map, Basemap, etc) in an ArcGIS Pro project using ArcPy.

However the ArcGISProject class does not seem to provide a method to do this.

I would like to be able to make sure that a particular Layout in my ArcGIS Pro project is the view that a copy of my project (made using the saveACopy method) opens to when I give it to someone.

Is it possible to do this by any other means from ArcPy?

If I were using ArcPy with ArcMap the equivalent would be accessible in the activeView property of the MapDocument class.


texts: ArcGIS Javascript API AMD require syntax and Dojo ContentPane compatibility

texts:
I am trying to recreate the ArcGIS Javascript API/Dojo ContentPane layout template shown here for my own custom application. I would like to switch from using dojo.require() legacy syntax to using AMD module require syntax.

Here is how I have structured my Javascript (removed everything but the basemap from the ArcGIS template):

I didn't change much about the HTML except for changing the periods to slashes consistent with the AMD syntax:

        This is the header section

subheader

                Content for the second tab

        this is the footer section

Making this switch has caused a very annoying problem in which the actual map typically only fills half of its ContentPane when the page initially loads. This problem can usually be fixed by resizing the window.

I'm aware that esri/Map objects are typically instantiated to the size of their parent container. I am not well versed in the nitty-gritty of asynchronous Javascript, but the fact that this is happening makes me wonder if the ContentPane and Map are loading independently of each other and the Map isn't getting the size info of the ContentPane.

Is this what's happening?

If so, how can I avoid having to resize the window to get the map to fill the ContentPane?


texts: How to create a circle vector layer with 12 sectors with Python/pyqgis?

texts: From a given coordinate and radius I need to create a circle with 12 sectors. The sectors should start from north (0 degrees) in 30 degree steps clockwise. I could also imagine to create the sectors by triangles, if that would be easier.


texts: How to change layer transparency without changing border lines?

texts:
I have a raster image of elevations. I have added a vector shapefile that contains local political divisions. I want to keep the elevation colors, but I want to keep the black lines representing the political divisions as well.

I have reduced the transparency of the political divisions shapefile, but I can't figure out how to have this NOT affect the political division lines. Advice for reducing the fill transparency while keeping the black border colors at full?


texts:
I'm trying to update the units used in labelling in QGIS 3.0.

The following answer shows how to do that but requires the user to overwrite all the changes:
Labelling line layer using QgsPalLayerSettings() QGIS failed

I want to edit only the label units and nothing else about the label style.

For example, I can edit the line size unit in the following way:

iface.activeLayer().renderer().symbol().setOutputUnit(2)

Is it possible to edit only the label size units (i.e MM/Pixel etc...) via Python?


texts:
I am trying to run ArcGIS 10.1 flow accumulation on a 10000 x 10000 cell unsigned 8Bit integer flow model (but 80% of that is no data). The data has been through a "fill" procedure once.

The only unusual aspect is that I have inverted the DEM, so the water flows from 0 into negative values.. could be an issue??

The model has been running for 24 hours, but while ArcGIS claims it is processing, the GPU shows no activity and the output file is empty.

My machine is a 2Ghz quadcore, 64 bit Win7 OS. How long might I expect this to run, and is there anyway to tell if it is working at all?

Thanks


texts: Removing lines that do not provide connectivity using FME

texts: I have two datasets one for lines and another for polygons, the lines data that contains lots of small line segments that are not needed. I would like to delete the lines only if they are not connected to another line or polygon and only leave in the lines that provides connectivity. I would like to do that using FME but I'm struggling to find the right transformers, any leads please?


texts: Calculating destination point using Python

texts:
I am trying to write a function that calculates the GPS (decimal) coordinates of a point that is far of a given distance from a given point.
Let say that we start from point 1 (The stripe of Dublin) and we want to move 25km east.

This is what I have written:

import math
# Coordinates of "The Spire" (Dublin)
lat1_dg=53.34972
lon1_dg=-6.26028
bearing=90 #degree from North
distance=25000 # distance in meter
alpha = math.radians(bearing)
radius=6371000 # radius of Earth (in meter)
gamma = distance / radius
lat1_r= math.radians(lat1_dg)
lon1_r= math.radians(lon1_dg)
Δlat_r = gamma * math.cos(alpha)
lat2_r = lat1_r + Δlat_r
Δlon_r = gamma * math.sin(alpha)/math.cos(lat1_r)
lon2_r = lon1_r + Δlon_r
lat2_dg=(math.degrees(lat2_r))
lon2_dg=(math.degrees(lon2_r))
print(lat2_dg,lon2_dg)

This outputs the following, which is close to the correct answer

  53.349720000000005 -5.883634969206713

but if I calculate the distance between point 1 and point 2 with this tool, I get 25.008 meters...why?


texts: Change Layer's Label Property "Display coded value description"

texts:
I created labels for layer in according with this example
They display coded value description from domain, but I need value itself. I can resolve it manually by unchecking this option

  "Layer's Properties -> Labels -> Expression... -> Display coded value
  description"

How can I do that programmatically with ArcObjects 10+?


texts: How to begin creating a linear referenced gdb with NHD

texts: Our legacy fish distribution data contains data by reaches. Reaches are arbitrary locations on streams, eg. reach 1, 2, 3. Each reach's from-to are known, but not by x,y coords. An example is 'brown trout are found in this stream reaches 1 and 2, not 3".  Our data has many fish species with such data. I want to create a linear referenced fish distribution database that can dynamically segment this data.  It seems I should create a route layer from the NHD for our reaches first, then create event layers for each spp. My question: "Is creating a route layer for reaches necessary? or can I digitize directly from the NHD?  I am new to LRS.


texts: Replacing NULL value with Zero in geodatabase table using Python parser of ArcGIS field calculator?

texts:
I am having some difficulties to change "NULL" values to zero in geodatabase. I have tried the following python code in ArcGIS 10.1 ("Field Calculator") to change it but did not work.

Pre-Logic Script Code:

def RemoveNULL(x):
    if x is None:
        return ''
    elif x == '':
        return '0'
    else: return x

Calculation field:

RemoveNULL(str( !DEAD_VOL_PER_HA_SPP1_125!))

Error message:

ERROR: 999999: Error executing function
The value type is incompatible with the field type. [dead_pine]
Failed to execute (calculate field)


texts: Carto legends are missing when embedded?

texts:
I've been using CartoDB at work to create maps and then embed them on our website. This has been working fine until recently, when the platform switched to Carto Builder.

Now, when embedding maps on the website, the maps are missing their legend although I created them and can see them on Carto. I also made sure to check that the legends box is ticked on the Map Options bar.

Why are the embedded maps missing the legend?


texts: NConverting Shapefile coordinates into lat/lon?

texts:
I'm having a little trouble importing a shapefile into PostgreSQL.  Using shp2pgsql, I'm able to import the non-geometry data just fine.  However, the geometry data is coming in based on some coordinate system rather than the lon/lat we require.  I'm hoping to make this conversion it the import stage, but I seem to be unable to find the magic combination of projection conversion parameters in shp2pgsql that accomplishes this.

Attempting to import using...

shp2pgsql -c -d -D -W LATIN1 -s 4269 -I 'path/to/shapefile.shp' my_table | psql -d mydb

...works except for the geometry column conversion I need.  Querying...

SELECT ST_AsText(geom) FROM my_table LIMIT 1;

...produces...

"MULTIPOLYGON(((-9619467.85573143 3856511.77685212,-9619466.74264784 3854145.2378255,-9619466.73474416 3854144.27900372 (...)"

My goal is to produce something like...

"MULTIPOLYGON(((-123.1 45.1,-123.2 45.2, -123.3 45.3 (...)"

(I'm using placeholder lon/lats, obviously).  The projection file for the shapefile reads:

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG",
"7030"]],AUTHORITY["EPSG","6326"]],PRIMEM
["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT
["degree",0.0174532925199433,AUTHORITY["EPSG",
"9122"]],AUTHORITY["EPSG","4326"]]

I'm at a loss as to what to try next.


texts:
I created a simple map with a point and a buffer (3 km) in QGIS2.18 (and also 3.4) with CRS WGS84 UTM33.

With QGIS Composer I created the layout adding the scalebars (for 500 m and 3 km) and printed a map in scale 1:50000.

After exporting the map as PDF and printed it, I checked the distance on the paper map  and I noticed that these are not so perfect: the measures of scalebars are ok (eg: 1 cm in the scale bar is equal to 1 cm on the ruler) but width and height of grids are little bit lower.

I also noticed that exporting the map as image (png) I have errors in scalebar (only in the 3 km one), buffer and grids, instead exporting in PDF, errors are visible only on the grids.

Pictures of measure are availabe at this link.

Did I wrong something in print setting or in the map?


texts: Create geometry from a result of intersection

texts:
Split a polygon in half with an extended line

public Ipolygon CreatePolyFromIntersection(Ipolgon polygonA, Ipolygon         bounradyThatcontainsPolylgonA, IFeature lmark){
            IConstructLine cLine = new LineClass();
            var parcelSegmentCol=(ISegmentCollection)PolygonA;
            for (int i = 0; i < parcelSegmentCol.SegmentCount; i++) {
            var crntSeg = parcelSegmentCol.Segment[i];
            crntSeg.SpatialReference = lmark.ShapeCopy.SpatialReference;
            cLine.ConstructExtended((ILine)crntSeg, esriSegmentExtension.esriExtendEmbedded);
            var line = (ILine)cLine;
            IPolyline extendedline = new PolylineClass();
            extendedline.FromPoint = line.FromPoint;
            extendedline.ToPoint = line.ToPoint;
            extendedline.SpatialReference = crntSeg.SpatialReference;
            ITopologicalOperator createInter = (ITopologicalOperator)extendedline;
            IElement em2 = new LineElementClass();
            em2.Geometry = extendedline;
            var pointsOnPolygonB= (IPointCollection)createInter.Intersect(bounradyThatcontainsPolylgonA, esriGeometryDimension.esriGeometry0Dimension);
            var point1=pointsOnPolygonB.Point[0];
            var point2=pointsOnPolygonB.Point[1];
            var intersectingSegment = new PolylineClass();
            intersectingSegment.FromPoint = point1;
            intersectingSegment.ToPoint = point1;
    // now I have points on Polygon B, all I have to do is is query those points on the polygon, get their distances on the polygon, and get the subCurve using those distances, I will get a curve, and with this curve I can add the intersectingSegment I created from those 2 points and I will have my polygon any comments ????


texts: Data conversion between Shapefile and Geopackage: does it affect accuracy of data?

texts:
I use Shapefile as interchange format and Geopackage as the format to to processing. So I convert file format at least twice between Shapefile and Geopackage.
Does/can this somehow affect the data accuracy and precision? I'm not asking about the limitations of Shapefiles like file size, field lengths etc., but accuracy of the data itself (comparable to reprojecting data several times: this can affect accuracy as transformations are not always exact). I work in QGIS, if that matters.
To resume a summary workflow:

Someone sends me data as a Shapefile
Than I convert it to a Geopackage to process it
Finally I convert it back again to a Shapefile to send it back

So to formulate the question otherwise: can I expect more accurate data if I avoid conversion and process the data as Shapefile?


texts: Changing date format in field name from d2000_5000 to Dddmmyyyy?

texts:
I am trying to change the name of a field by using Python code where the current name is in the format (e.g.) d2000_5000, and 5000 is the decimal year (0,5 y) added to 2000. So i think that this can be done with the Alter Field tool, by extracting characters, then doing 365*5000/10000 to obtain days, and then adding these days to the initial date of 2000: 2000-01-01 + days (within a date sum that takes into account years with 29 days in Feb).

Somebody can help me proposing a code snip to do it?


texts: Reprojecting a raster image using gdalwarp

texts:
I am trying to reproject an ECW Orthophoto from 25830 to 23030 srid.
I tried the following but I am getting an error:

C:\Program Files (x86)\FWTools2.4.7>gdalwarp -s_srs 25830 -t_srs 23030 C:\orto.e
cw C:\orto_ed50.ecw
ERROR 1: Translating source or target SRS failed:
25830


texts: Can I use google maps navigation within my app

texts:
I am considering building an app which includes google maps. The app functions by tracking the GPS location of the vehicle and showing it to the respective user on the Android app with some additional stats. I also would like to include turn by turn navigation. Based on what I have read in Google's terms and conditions it does not allow to use navigation but I have seen on the app runtastic on google play store using google maps navigation. I wanted to know how can I avail this functionality. Can it be made available to use on purchasing of the work API? Also does Google provide it's SDK for the same.

It's been a few days now that I had sent this inquiry to google but haven't received any response so decided to post it here.

Source: https://developers.google.com/maps/terms#section_10_4 / c. / iii.

  10.4 Restrictions on Unfair Exploitation of the Service and Content.

  ...

  c. Restrictions on your Maps API Implementations.

  ...

  iii. No navigation. You will not use the Service or Content for or in connection with (a) real-time navigation or route guidance; or (b) automatic or autonomous vehicle control.


texts:
I've just installed Geoserver 2.3-RC1 on windows 7 64bit and JRE 6.43 and downloaded the relevant sql extensions and Microsoft files .
As instructed , I've copied "gt-jdbc-sqlserver-9.0-RC1.jar" to "C:\Program Files (x86)\GeoServer 2.3-RC1\lib" and also the x86 versions of "sqljdbc_auth.dll" and "sqljdbc_xa.dll" to "C:\Windows\System32".

I've restarted geoserver and went to Add Store and still didn't see SQL option.

I've tried doing the same with the x64 versions of the dlls , same result.

What am I missing ?!


texts: Error on domReady callback: TypeError: b is null arcgis

texts:
I'm getting this error "Error on domReady callback: TypeError: b is null" when including arcgis javascript api to my project.

<title></title>

This script is in page called map.aspx and I request it using iframe like this

I'm not sure what is the problem, any help?


texts: Can I compare Moran's I z-scores across areas?

texts: I would like to compare the spatial clustering of elderly people (measured as percent of population age>64 in each census tract) across five cities. I have run "Spatial Stat. tools -> Spatial autocorrelation") and obtained a Moran's I z-score for each city. I have applied exactly the same tool parameters for the five cities. My question: does it make sense to compare Moran's I z-scores for the five studied cities? Or can I say that the city with the highest z-score has the greatest level of clustering?


texts:
I have been using the shortest path (Point to Layer) option but sometimes I get "There is no route from start point (x, y) to end point (z,w)".
After checking my network, I notice that this message happens usually if there is a smaller network not connected to the largest network that is closer to either the start or ending point.  Is there a way to prevent this from happening?
Example in the image. Area 96 has a shortest path found by the algorithm shown in yellow, which is around 300 meters from the ending point.
For area 97 there is no shortest path found.  The area in blue is a small network without connection to the origin point.  The area in orange is connected to the larger network that connects to the origin point.  THe blue area is closer to the 97 point than the orange point, however the "orange" line is 200 meters away from point 97, which is closer than the distance that gave me a correct shortest path for zone 96.


texts: Why won't my street addresses georeference?

texts:
I've been using Carto for about a month or two now and it had been working perfectly, however it suddenly malfunctioned today. I updated my excel document (nothing changed except the number of rows), uploaded it to Carto, and then analysed it for Street Addresses. For some reason the street addresses won't analyze and when I try to style the layer it says "there's no geometry in your data that could be styled".

The thing is I've updated this document about 30 times before, uploaded, analysed, and styled it and it worked every time. I had uploaded, analysed, and styled a layer about 5-10 minutes before this, so it was working fine until about 1:00pm EST. Is there anything in Carto I could have clicked that would stop it from reading the addresses correctly? The excel document is still in the same format as before, and I've uploaded files equally as large, if not larger, before so I know it's not a size thing.


texts: QGIS creating a layer to draw over Google Streets Map

texts:
Is there a way to create a layer that covers Google Maps blocks

for example this is the google map,

and this is the layer that i want to create

since i cant use  snap on google maps im trying to work around it and find a way to automatically create those blocks


texts:
I know this question has been asked a million times already, but I can't find many detailed responses on exactly what I want to do.

As far as I can tell, the full C code for the fuller projection is available here: http://www.rwgrayprojects.com/rbfnotes/maps/graymap6.html

Since I'm not using it for any commercial purposes (personal use only) it doesn't violate the licensing. It would appear as though it would only need a few changes (if any) to incorporate into proj4. My only question is, how do I incorporate it into proj4 and have it recognized by gdal? Is it as simple as putting the ported, proj4 compatible version into a PJ_Fuller.c file and building the project? How would I assign it my own custom ESPG code, and once proj4 is rebuilt, would I just symlink my custom branch of it so that QGIS looks there and sees it? Or would I also have to rewrite the proj4 file in c++ for gdal?

I'm running Debian if that makes any difference.


texts: Support ArcGIS Enterprise + Portal across multiple domains

texts:
I have an ArcGIS Enterprise + Portal environment where I would like to host content across two domains to two separate sets of users (domain1.com and domain2.com).  I would like to find any documentation (case studies, architectures, best practices, etc) that would help me architect a solution that minimizes data redundancy and prioritizes security.  Both systems are on Windows server using Active Directory (not Azure AD).

I would like to use Active Directory as the identity store.

Does this mean I need two environments for each domain or is there a recommended and tested way to support two domains on one cluster?