texts: Calculating spatial autocorrelation of raster values?

texts:
R is perfectly willing to compute a Moran's I statistic (and probably all of the other autocorrelation metrics) for a raster.  As best I can tell, PySal and ArcMap's autocorrelation tools specifically require feature data, but it's easy enough to game that requirement by converting a raster into a set of regularly spaced points or a random point sample.

My question is: should I?  And if not, what to do instead?

Specifically in reference to this earlier answer:

  I would point out (pardon the pun) that ecological phenomena, such as
  vegetation cover or habitat type, do not meet the assumption of a
  point process. A point pattern (process) represents a set of "events"
  that exhibit random or spatial characteristics that can be quantified
  statistically. The idea that a simulated point process is a sampled
  representation of a empirical process occurring across discrete areas
  is an erroneous assertion. Discrete areas do not represent a point
  process and point pattern analysis statistics are quite inappropriate!

If rasters are the wrong material for the point pattern autocorrelation measures, what is the correct way to quantify the degree of spatial autocorrelation -- local and global -- in the raster's values?

Elaboration:

This work concerns the statistical effects of raster smoothing/aggregation.  The original rasters contain either elevations or forest canopy heights taken from lidar, distance from a linear feature, counts of conifers, or (still looking for a good example of) categorical data, as well as some wholly artificial rasters designed to embody various spatial structures -- random, ordered transition from high-to-low, isotropy, etc. These are then aggregated by different statistics to a range of cell sizes and shapes.  Block statistics so far, with an eye on focal statistics later on.  I would also anticipate volume or biomass as a future value, if that requires something different.


texts: How to define the angle of points on slope?

texts:
I have point data set and DEM and my task is to identify how many points are falling on the slope of more than 40 degrees.
How could I do this in ArcMap?


texts: How to install two versions of QGIS on Ubuntu 12.04?

texts: I would like to run multiple versions of QGIS on Ubuntu 12.04. Can I do that, and how if so?


texts: Identify optimal subset of polygons based on cumulative overlap with a separate layer?

texts:
Take the example below, where there are four shapes overlapped by a grid of points. Suppose that each shape has a "score" attribute linked to it, where score = number of sides of the shape (so circle = 1, triangle = 3, etc.). I want to find the optimal set (S) of k points (where k is a number decided beforehand), where S, when treated as a single unit, overlaps polygons with the highest possible summed score.
Provided k = 2 in the example below, the optimal set S would contain points 3 and 5, which collectively overlap shapes with a total score of 17 (hexagon [six] + rectangle [four] + square [four] + triangle [three]), higher than any other possible pair of points. Note that the answer is not points 2 and 3, because the score from the hexagon cannot be counted twice (so the total score of points 2 and 3 would be 13 [6+3+4] rather than 19 [6+3+6+4]). Thus, scores can't be calculated for each point individually using spatial join; they can only be calculated for specific pairs of points based on the unique shapes that the pair overlaps with.
Is there a way to have ArcMap determine the optimal set of points (S)? I also realize that this can be figured out manually for simple examples, but in my case there are thousands of points, which overlap range shapefiles for thousands of species, with each species having an associated score. My goal, specifically, is to identify the 20 points on the map that have the highest cumulative score-weighted species richness possible. I can't just choose the 20 points that have the highest scores individually, because they're all clustered next to each other on the map and thus harbor almost exactly the same species, which would limit the cumulative biodiversity captured by the top 20 points in total if using this method.


texts: Selecting and clipping layers within user defined area using PyQGIS

texts:
I am new to QGIS and am trying to write a script that will allow the user to define a region on the map (either polygon or freehand) and clip all displayed map layers within this region to the user-defined selection area so that line lengths can be calculated and totaled within the region for each layer.
Using the example below, suppose there are two vector layers of polylines. The long vertical line is on one layer, and the pink lines are on the 2nd layer.
Instead of selecting the entire intersecting lines from endpoint to endpoint, I need to select/clip only the part of the lines that are within the polygon boundary.

The output I am looking for would be something like:
Layer Name 1: Total Lines Length: 55
Layer Name 2: Total Lines Length: 145
Layer Name 3: etc.

I would store this output in a table or csv file.
I am not sure how to get the user defined polygon region as input, how to get the displayed layers that fall within that boundary, and how to clip their extents to the region so that I can get their clipped lengths.
d = QgsDistanceArea()
d.setEllipsoid('WGS84')

#this is just selecting the active layer,
#but I need to get layers from boundary

layer = qgis.utils.iface.activeLayer()

"""
what I need is something like this <pseudocode>:

selection_layers = interface.CurrentSelectionPolygon()
for layer in selection_layers:
    new_temp_layer = clip_to_selection(layer)
    for feature in new_temp_layer.selectedFeatures():
       line_length = d.measureLength(feature.geometry())
       #add line length to total for layer etc.

"""

#getting the user selected features from active layer just to test
features = layer.selectedFeatures()

#
lengths = []
if len(features) != 0:
    for f in features:
        geom = f.geometry()
        #use measureLength to get accurate projected values
        length = d.measureLength(geom)
        #convert this to feet
        c_length = d.convertLengthMeasurement(length, 2)
        #0 is meters, 2 is feet
        lengths.append(c_length)

# this is outputting the sum of the single layer's line lengths
total_length = sum(lengths)
print(f"total length: {total_length}")

I cannot figure out how to get the user's map drawn selection polygon to use as a clipping mask.


texts: I need to interpolate lat/lon points along a small circle, with a given centre and radius.  I've successfully used these formulas for great circles.  What is the equivalent for small circles?


texts: How to query osm2pgsql tables for address geocoding?

texts:
I am new to geocoding and am trying to build a geo application. I used osm2pgsql to convert my osm data to pgsql gis extended database with the following tables planet_osm_line, planet_osm_polygon, planet_osm_point, planet_osm_ways, ...etc

My first query is to find the distance between 2 randomly given addresses using the above generated tables.

Can someone explain how to query these tables to attain the result ?

The duplicate Convert OSM File to get addresses using Osm2pgsql which was linked does not have proper information to query tables.

For example, I need to know which column has the name of the streets and which column has the house numbers and which column has names of the city or state and how to query between them to geocode or reverse geocode.


texts: Determining if specific lon/lat in various coordinate systems?

texts:
I'm a newbie in Gis matter and don't want to make it very complex.
I want to plot a point somewhere and know it's lon/lat in some coordinate systems indicated by EPSG code.

I'm wanting to know if there is a website that provides a map, that we can plot a point on, and the website tells me for example: this point is (lon=123456.123456, lat=456789.456789) in EPSG 32639 AND (lan=147.147, lat=369.369) in EPSG 3857 AND ... ?

Or at least tells me What EPSG is mostly used for this area/country?


texts: Labels are not showing in a published ArcGIS Feature Service?

texts:
I'm using ArcMap 10.1 to publish a feature service.  The labels show on the mxd, but not on the service (on a web application or on mobile application).

Do I have to create an annotation layer or am I just not setting something up correctly in my mxd?


texts: Does GDAL support multi-process ReadAsArray

texts: I have a very large TIFF file that I need to read and process block by block through ReadAsArray. If I have multiple processes calling this function at the same time, can I return the correct result?


texts: List with Strings as Filter for arcpy.mapping.ListLayers(mxd, "", df)

texts:
In my MXD with 14 different Raster-Sets (12 Months(January, February...), 1 GRID, 1 Border-Raster) I am looking to create a list (Layer Objects) with the function

arcpy.mapping.ListLayers(mxd,"FILTER",df)[0].

I don't know how to specify a filter (Wildcard) to create a list containing only the 12 Months. So far I tried to filter the list monlist with a tmplist containing all the months as Strings.
In the following example, the line

monlist = arcpy.mapping.ListLayers(mxd, tmplist ,df)

does not work. In result I get an empty list monlist=[]

import arcpy
raster = arcpy.GetParameterAsText(0)
outPDF = arcpy.GetParameterAsText(1)
hill_bool = arcpy.GetParameterAsText(2)

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
hillshade = arcpy.mapping.ListLayers(mxd, "hillshade", df)[0]
hillshade.visible = hill_bool

tmplist = ["Januar","Februar","Maerz","April","April",
"Mai","Juni","Juli","August","September","Oktober","November","Dezember"]
monlist = arcpy.mapping.ListLayers(mxd, tmplist ,df) #Select only the months
addLayer = arcpy.mapping.Layer(r"C:Python\Legende_2.lyr")
arcpy.mapping.MoveLayer(df, monlist[:-1],hillshade, "AFTER")

for element in monlist:
    element.visible = "TRUE"
    arcpy.ApplySymbologyFromLayer_management("","addLayer")
    if hillshade.visible = "TRUE":     #<<< ???
        element.transparency = 35
arcpy.RefreshActiveView()
...

#...align Legend and other objects -> ExportToPDF
arcpy.mapping.ExportToPDF(mxd, outPDF)


texts: Join a grid to points that are close but aren't overlapping

texts:
I have a point grid that overlaps points on a map. I want to join this grid to the points in order to keep only the parts of the grid that contain a point.

However, there are some points that lie between the circles of the point grid (in the empty space) and therefore aren't counted as overlapping and aren't included in the output.

The whole point of my grid though is that it should cover the entire area. I know I could get it to join with these points by making the grid a square grid (and eliminating the empty space) or by making the circles of the point grid so large that they fill most of the empty space. However, I want to keep the circles the same size.

Is there a way of joining the circles of the point grid to a point that is NEARBY rather than actually touching?


texts: Saving MapCanvas with RubberBand items as image

texts:
I want to save map canvas as image (Project > Import/Export > Export Map as Image...) together with rubber band items(colored areas in the image). But QGIS just saves features of shown layers on map canvas.

I get the same result when using PyQGIS.

iface.mapCanvas().saveAsImage(filename)

How can I save QgsRubberBand items by PyQGIS?


texts: Mapping Geocoded addresses onto Census shapefiles in R

texts:
I have a list of addresses that I have geocoded to get the corresponding longitude and latitude.

> AddressDetail <- geocode(match_address)
> coordinates(AddressDetail) <- c("lon", "lat")

Now I believe I should assign a projection so that R knows how the lon/lat coords are set up

> proj4string(AddressDetail) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

Now I read in the shapefile that uses the NZMG projection

> Shape_MB = readOGR("/Work in progress/MB06_LV2/", "MB06_LV2")

Now this is where it starts to go wrong. I try transforming the projection of the lon/lat addresses into NZMG projection.

> AddressDetail <- spTransform(AddressDetail_v2, CRS("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs")
non finite transformation detected:
lon lat
Error in spTransform(xSP, CRSobj, ...) : failure in points
In addition: Warning message:
In spTransform(xSP, CRSobj, ...) : 9 projected point(s) not finite

So when I continue on to use the over function to find the corresponding census meshblock index of the polygon that the address is within I get the following error.

> temp <- over(AddressDetail_v2,Shape_MB,MB06)
Error: identicalCRS(x, y) is not TRUE

I'm a bit of a novice with geospatial stuff in R so would appreciate any advice!


texts:
I'm quite new in ArcPy.

I'm able to get a list of values from a Layer Feature and for specific field, however I would like to edit that attribute after is displayed on the screen. Because the value is inside the list if I type another value is telling me is not in range in the list, so raises an error.

Is it possible after geting the list of values modify directly into the parameter value?

I was thinking that if the value is not in the list, then type your new value as an input on the same parameter box without raising an error.

This is my Validation Toolbox python code at the moment:

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""

  def updateParameters(self):
    fc = "my_feature_class"
    rows = arcpy.SearchCursor(fc)
    if self.params[0].value:
      self.params[0].filter.list = sorted (list(set(r.getValue('my_field')for r in rows)))
      del rows
    else:
      self.params[0].filter.list = []
    return


texts: Comparing geometries with having a buffer

texts:
I have a small road portion of a country in a form of linestrings. I want to compare it with the whole area of country(also linestrings) and get the features of the matching linestrings. The matching should be so exact that no extra points from whole country should appear in matching linestrings that does not appear to be in the small road portion. I have been comparing them with ST_intersect but the area I got contains some extra linestrings from the whole Area of country that are not in small portion of road. The techniques I have been reading suggests that I should be creating a some kind of a buffer to create the comparison. Somebody know an technique to get the exact comparison of these two linestring?. I have just started using POSTGIS so I am newbie.

The query I am using is this one but contains extra lines from WholeAreaTable which are not in SmallPortionTable

CREATE Table ExtraLineArea AS select id ,link.Geometry from WholeAreaTable WAT, SmallPortionTable SPT where
ST_Intersects(WAT.Geometry,SPT.Geometry)

I am selecting the id and geometry of WholeAreaTable when both geometries intersect. In the image brown line represents that intersecting area of both geometries. As WholeAreaTable contains all the geometries of country it is giving some extra line other then intersecting area which is purple line in the image(see on top of image single purple line).

I want comparison of WholeAreaTable and SmallPortionTable in such a way that this purple line should not appear in the final result.


texts: Can distance be calculated to a destination feature based on an attribute in the source?

texts:
All,

Is it possible to specify which destination point to calculate distance to based on an attribute in the source?

STDistance will give me the distance each student (point) is to their closest school (point).

When I pull a list of students based on their boundary school number (the school they should go to based on their proximity to their house), I also get back their choice school number (the school they actually go to).

I would like to calculate the distance each student is from their choice school.

I'm using SQL Server 2012, hence the STDistance function (rather than ST_Distance in PostGIS)

Any thoughts appreciated!

-m


texts: Opening images georeferenced (.aux.xml) in ArcMap in QGIS?

texts:
We have several ArcMap georeferenced images (raster datasets) that we are attempting to open in QGIS without any luck.  These are not inside of a geodatabase and are not  Arcinfo Binary Grid format or Arcinfo Binary Coverages.  QGIS supports GDAL rasters and the ArcMap created .aux.xml file contains the georeference points, but QGIS apparently does not recognize this file as it places my Texas image in a very small area of the Atlantic ocean off the coast of Africa (with a google maps basemap added).  (THe only other files in the raster dataset are .tfwx and .ovr)

Is there something special I need to do within QGIS to apply the coordinates in the aux.xml created by ArcMap, or something else I'm missing?

Here are the contents of one such aux.xml file:

<?xml version="1.0"?>

-<PAMDataset>

<Metadata/>

-<Metadata domain="Esri">

<MDI key="PyramidResamplingType">NEAREST</MDI>

</Metadata>

-<Metadata domain="xml:ESRI" format="xml">

-<GeodataXform xmlns:typens="http://www.esri.com/schemas/ArcGIS/10.1" xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:type="typens:PolynomialXform">

<PolynomialOrder>1</PolynomialOrder>

-<SpatialReference xsi:type="typens:ProjectedCoordinateSystem">

<WKT>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]]</WKT>

<XOrigin>-22041259.177065704</XOrigin>

<YOrigin>-33265069.154224794</YOrigin>

<XYScale>135385247.70505452</XYScale>

<ZOrigin>-100000</ZOrigin>

<ZScale>10000</ZScale>

<MOrigin>-100000</MOrigin>

<MScale>10000</MScale>

<XYTolerance>0.001</XYTolerance>

<ZTolerance>0.001</ZTolerance>

<MTolerance>0.001</MTolerance>

<HighPrecision>true</HighPrecision>

<WKID>102100</WKID>

<LatestWKID>3857</LatestWKID>

</SpatialReference>

-<SourceGCPs xsi:type="typens:ArrayOfDouble">

<Double>5.994281757832141</Double>

<Double>5.6742238617780458</Double>

<Double>2.0592992159490677</Double>

<Double>4.5229583592636118</Double>

<Double>2.014907522695573</Double>

<Double>8.8668953443025202</Double>

<Double>2.8272115205686532</Double>

<Double>7.7051733542912189</Double>

</SourceGCPs>

-<TargetGCPs xsi:type="typens:ArrayOfDouble">

<Double>-10752195.345278211</Double>

<Double>3857531.2694506948</Double>

<Double>-10756947.906142347</Double>

<Double>3856205.4206410511</Double>

<Double>-10756935.616529588</Double>

<Double>3861470.3939956115</Double>

<Double>-10755980.997953685</Double>

<Double>3860058.8390891682</Double>

</TargetGCPs>

<Name/>

</GeodataXform>

</Metadata>

</PAMDataset>


texts: Near tool producing odd results for batch-snapping; how to best perform batch-snapping?

texts: I am trying to use the Near tool in Arc10 as described in Snapping points to a line using ArcGIS geoprocessing tools for the purposes of batch-snapping a feature class of points to a feature class of lines. The XY coordinates produced by the Near tool do not exactly map onto the lines - there can often be about a foot distance between the line and the new point. The new points also do not appear to be on the shortest-distance path from the original points. The projections are the same for the two feature classes and the data frame, but points are in meters while lines and data frame are in feet. In attempting to fix this problem, I have projected my original points into feet before using Near, with no change in the outcome. Any suggestions on why this method doesn't work for me, and how to get my points to batch-snap to lines, are greatly appreciated!


texts: After the qgis2web plugin was upgraded last week the popup window shows all fields in the attribute table when you click on a feature regardless of which ones were selected in the dropdowns when creating the webmap.  Is it just me or is anyone else experiencing this?


texts:
First of all, I want to apologize because I am a beginner in this area and not a native English speaker, so sorry if I ask simple questions or if I am hard to understand.

So I have DEM files in the form of ASC files representing squares of 5 km and a resolution of 5 m (so 1 million points/file). What I want to do is import them in Blender to edit them and study them. I am able to convert the ASC files to TIF files and import them to as points cloud with the BlenderGIS plugin but then I need to skin them with a plugin which make mistakes and the files are far too large to be corrected manually. The other method i tried was to import the file in Blender "As DEM" and it work great but when I switch to edit mode, the mesh is flat and I cannot edit it the way I need. Do you see any method or tool that I can use to help me ?

As I said, I am a real beginner and I may have missed something obvious so if you think of anything please try to help me !


texts:
I am using GEE and I'm pretty new with this engine. I'm trying to download several images from Landsat collection. I provide the script for you.
Here is my question:
How can I download each NDWI image as a separate file (rename it with the date and path&row) and export them into my Google Drive so I can use them with different software?

Map.centerObject(geometry2);
var landsat = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate('2019-03-01','2019-04-30')
.filterBounds(geometry2);
print(landsat);
Map.addLayer(landsat,{bands : ['B6','B5','B3'],min: 0, max: 1});

var NDWI = landsat.map(function(img){
  var index = img.normalizedDifference(['B3','B5']);
  var clip = index.clip(geometry2);
  var id = img.id();
  return clip.rename(id);
});

print(NDWI,'NDWI');

var stack = NDWI.toBands();
print(stack,'stack');
Export.image.toDrive({
  image : stack,
  description : 'NDWIstack',
  scale : 30,
  region : geometry2,
  maxPixels : 1e10});

This script maps the NDWI function on all images in the specific location and within the specific time range (from March 2019 to the end of Apr 2019). The result is 395 NDWI images. When I stack them (using toBands() ) all images become ONE image with 395 bands (am I right?) When I export them I find the images with wrong names (please note that I am using the ID of each image to rename the output images). The images are also more than the results (i.e. 395). For example for a very smaller region (geometry) the results is 9 images. When I export them into my Google Drive the images are 42! and of course with names like "NDWIstack-0000011008-0000055040.tif". I downloaded one image and I cannot open it in any image processing software (ArcMap, ENVI)


texts: Create 10.0 compatible Mosaic Dataset with 10.2 (or something else)

texts:
I have a client who requested a ArcGIS mosaic dataset, but as I found out today they only have licenses for ArcGIS 10.0, and I only have licenses for 10.1 and 10.2. The mosaic datasets I create cannot be viewed on their systems.

Is there any way to create a backwards compatible mosaic dataset with what I have? If not, are there any tools available elsewhere that would accomplish this (aside from getting a license for 10.0)?


texts: Returning US Census geographies within range of coordinate using PostGIS and Tiger Decoder?

texts:
TL;DR:

How do I write a query in SQL (postgres, postGIS, tiger extension) that will return census geographies (blocks, tracts, zcta) within an X mile radius of a lat/lng coordinate?

Background:

I would like to feed a query a lat/lng pair and receive back census geographies (block, blockgroup, tract, county, state).

EDIT
Piggybacking off this answer, I have a query that will return the information I need at a single coordinate. How do I extend it to include a radius?

Finding census block for given address using Tiger geocoder

SELECT tabblock_id
FROM tabblock
WHERE ST_Contains(the_geom, ST_SetSRID(ST_Point(-71.101375, 42.31376), 4269))

This answer here suggests I can use something like this to get the data I need. With the query below and the previous one, I think can mash them up to get the radius information I need.

Return all results within a 30km radius of a specific lat/long point?

SELECT *
FROM your_table
WHERE ST_Distance_Sphere(the_geom, ST_MakePoint(your_lon,your_lat)) <= radius_mi * 1609.34

However, now I'm seeing my tabblock table is empty.

geocoder=# select count(*) from tiger.tabblock;
 count
-------
     0
(1 row)


texts: Removing noise from point cloud with lidR package

texts: I have unclassified Terrestrial Laser Scanning (TLS) point cloud data which has some noise. Is there a function in the lidR package which can be used to remove noise?


texts: Unexpected problem with some CSS styles on GeoServer for OSM data

texts:
My layerset went wrong after mapproxy caching. Some layers are not rendering, and one of them shows this error message when I want to watch a preview:

<ServiceException>The requested Style can not be used with this layer. The style specifies an attribute of water and the layer is: test_osm:water </ServiceException>

</ServiceExceptionReport>

CSS code of this style:

*{
  fill: #91f1ff;

}

[@scale > 5000001] [name='Байкал' or water='sea' or natural = 'coastline']{
 fill: #91f1ff;

}

[@scale < 5000001] [@scale > 500001][water='sea' or water='canal' or water='reservoir' or water='canal' or water='lake' or water='cove' or water='river' or natural = 'water' or waterway ='riverbank']{
 fill: #91f1ff;

}

[@scale < 500001][water='sea' or water='canal' or water='reservoir' or water='canal' or water='lake' or water='cove' or water='river' or natural = 'water' or waterway ='riverbank' or water = 'pond']{
 fill: #91f1ff;
}

What could be the problem? Validation is OK and CSS files really exist on workspaces/myworkspace/styles...


texts: How to run GDAL vector (OGR) tools from inside QGIS?

texts:
I want fingertip access to the GDAL/OGR tools, particularly ogrinfo, from the main QGIS screen.

I currently run ogrinfo.exe in a DOS shell but this seems clumsy.

Installing the GDAL Tools plugin only accesses the raster tools. (The User Guide confirms this.)

But in my Windows QGIS folder there is a file:

 C:\Program Files\QGIS283\apps\qgis-ltr\python\plugins\processing\algs\gdal\ogrinfo.py

Does this perhaps somehow let me run ogrinfo from the Python console?


texts: EPSG:4038 Gridset does not work with caching on GeoServer 2.2

texts:
I am trying to setup a cache for a GeoTIFF base layer. I set it up and I went to the built in GeoWebCache page where it generates test map for the caching. I saw the built in projection
- EPSG:900913
- EPSG:4326

I tried these two links and it shows open layers map and it generates cached tiles properly. So I set to work to create a gridset for ESPG:4038 using this manual: http://docs.geoserver.org/latest/en/user/webadmin/tilecache/gridsets.html#webadmin-tilecaching-gridsets (I added 20 levels of zoom)

I added this gridset to the base map layer cache and now I see these three items in the GeoWebCache page.

EPSG:900913
EPSG:4326
ESPG:4038

I tried the first two links and it works properly (showing the maps and I see the indicators that certain amount of disk space are being used by the cache layer). But when I try to open ESPG:4038, the generated OpenLayers preview map is blank and there's no activity in the disk space indicator for the cache. When I check on Google Chrome, it return 256x256 transparent png tiles.

I'd appreciate if anyone can point out which part I do wrong


texts: QGIS 3.4 Rule Based Symbology - Exact Matches Integers in String e.g. 2,12 + 1,2,3,12

texts:
I have a shapefile containing a String Type Field called 'Class'. This field contains one or more integers separated by commas.

I want to create Rule-Based symbology matching only 2, 12, 1 respectively.

I've tried the following rules with no luck, for the string '1,2,12' each rule matches both 2 + 12 e.g.

Rule 1 - 'Class' LIKE %2%

Rule 2 - regexp_match('Class','2')

I was thinking adding * before and after each integer when the field is being populated but that would make the field very long. Was hoping there would be a fix for this problem.


texts:
I started studying the examples provided in the "QGIS Python Programming cookbook".

I downloaded the point vector layer as example from https://code.google.com/archive/p/geospatialpython/downloads?page=2 and uploaded in QGIS.

The code below should append a point to the existing points of the vector layer:

vectorLyr = QgsVectorLayer('/qgis_data/nyc/NYC_MUSEUMS_GEO.shp', 'Museums', "ogr")
vpr = vectorLyr.dataProvider()
pnt = QgsGeometry.fromPoint(QgsPoint(-74.80,40.549))
f = QgsFeature()
f.setGeometry(pnt)
vpr.addFeatures([f])
vectorLyr.updateExtents()

But, at the end, any new points has been added.
May you help me?


texts:
I have a DEM raster with some small and large areas of NoData. I also have a raster defining zones (the zones consists of a rectangular grid).

I want to create a new raster marking the zones where there exists any nodata values in the DEM raster.

My initial approach was to use the Zonal Statistics tool, using any statistics type (for example Maximum) and uncheck the "Ignore NoData" option. According to what I can read in the documentation, this should give me an output raster with NoData in the entire zones where the input raster contained any NoData pixels. (I could then perform some trivial operations to give me the desired final result).

However, this doesn't work as expected. The output raster from the Zonal Statistics tool have some holes, but there are several places where there are definitely NoData in the input raster and data in the output.

So, are there any other possibilities of achieving what I want? Or have I misunderstood how I should use the Zonal statistics tool?


texts: Importing numeric fields from CSV to DBF in ArcGIS Desktop?

texts:
I am trying to import a CSV into a DBF using ArcCatalog (10.1).
In the field map, pictured below, I can only set an ID field to Long, Text, Double, etc. But there is no Numeric type as in DBF and shapefile attributes.

If I import the ID field as a Long, values will be truncated. If I import it as Text, it will be different type than Numeric. In either case, I won't be able to join the imported DBF to a shapefile with the ID.

How should I choose the type so that the ID in the imported DBF can be joined with a 'Numeric' ID in shapefiles?


texts: Marker icons overlapping street names in Tilemill?

texts:
I'm new to CartoCSS and I'm trying to properly add subway markers on a map using Tilemill and OSM data.
My problem is that some street names stay behind my marker icons :

I have add "text-allow-overlap: false", "text-placement-type: simple" and "text-label-position-tolerance: 10" on my road labels in order to make them avoid my markers but I still get these superimpositions.

This is my road labels code:

#mainroad_label[type='primary'][zoom>12],
#mainroad_label[type='secondary'][zoom>13],
#mainroad_label[type='tertiary'][zoom>13] {
    text-name:'[name]';
    text-face-name:@sans;
    text-placement:line;
    text-fill:@road_text;
    text-halo-fill:@road_halo;
    text-halo-radius:1;
    text-min-distance:60;
    text-size:11;
    text-allow-overlap: false;
    text-placement-type: simple;
    text-label-position-tolerance: 10;
    [zoom>=19] { text-size:12; }
    [zoom>=20] { text-size:14; }
}

When I set "marker-allow-overlap" to "false" it hides the markers located on text (street or place names). However, I'm wondering why it is not working in the other way ?

Here is my markers code:

#stations[Description='m'][zoom>=14] {
    marker-file: url(/Users/alexis/Documents/MapBox/Project/OSMBright/Metro-M.svg);
    marker-height: 11;
    marker-allow-overlap: true;
    marker-ignore-placement:false; }

Is it possible to make text disappear when located on a marker ?


texts: How to check the encoding in the .dbf files?

texts:
I am having difficulties when I try to load a shapefile into my database, the is because I do not know the encoding of the dbf portion of the shapefile.
The attribute tables are in Spanish, I have tried UTF8 and LATIN1/ISO-8859-1.

How can I check the encoding of the dbf file? I do not want to waste anymore time trying all the encodings.

The database is Postgres 9.4.4 with Postgis 2.2.1, I use Qgis 2.10.1 to visualize the data before loading it to the database.

Looks like the files were made in a mac, not sure though.


texts:
I have a .tif file (about 380 mb) that I am trying to load in R and plot.

file<-raster('tiffname.tif')

It is a formal class Raster Layer with the following information

class       : RasterLayer
dimensions  : 40000, 40000, 1.6e+09  (nrow, ncol, ncell)
resolution  : 0.00025, 0.00025  (x, y)
extent      : 80, 90, 30, 40  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84+towgs84=0,0,0
data source : \Data\trial\tiffname.tif
names       : X40N_080E_mtc02

Why are there no values for this tif file? Is that why I am unable to plot it?


texts:
Can anyone spot what is wrong here? I've selected single feature, the existing rings are obvious, it's not intersecting anything (even a feature in another layer...

http://johnsankey.ca/ringerror.jpg


texts: Define multiple fields to copy in graphical modeler in QGIS 3.6

texts:
I am looking to find the correct syntax for copying a defined set of fields when using the "Join attributes by field value tool" in the graphical modeler in QGIS 3.6.

I would like to know how the syntax for this when selecting multiple fields as I cannot seem to work it out. For example if I only wanted to join "Field_1" and "Field_2".

Please see my image below:


texts: Converting QGIS styles to ArcMap format (.style)

texts:
Background: Most of the people in my company work in QGIS, however there are some deliverables such as SSDM - Seabed Survey Data Model that have to be produced in ArcMap. As a part of SSDM we're supposed to supply the symbology, in .style format.
Problem: Manually creating each symbol used in QGIS in ArcMap is a pain, so I'm wondering whether there is a way to convert/translate the QGIS styles to ArcMap format? Preferably straight into .style file, but even loading each symbol one by one would be much faster.
I've done some research, I know that there is a SLYR plugin and several other solutions, but they all seem to be only converting Arc styles to Q. I haven't found anything about doing it the other way around. Any help will be appreciated.


texts: Difference between mask and bbox in geopandas read_file

texts:
Geopandas read has as options bbox and mask.
https://geopandas.org/reference/geopandas.read_file.html
What is the difference between those two options?


texts: How does one convert a .kmz or .kml to a shapefile using arcpy? tried the KMZ to Layer but it outputs a .lyr and a .gdb containing one feature class. I need a shapefile, or a way to convert the contents of the .gdb to a shapefile.


texts: Is there a straight way to copy or clone a QgsVectorLayer Object?

texts: Is there a straight way to clone or "deepcopy" a QgsVectorLayer Object with python? Unfortunately the python module "copy" is not successful with that.


texts: Viewshed plugin results in memory error in QGIS 2.6

texts:
Everytime I use the viewsheds plugin up to this point I get an memory error. I projected my raster files and the point which I added as an observervation point into a metric system (EPSG:3857 - WGS 84 / Pseudo Mercator) and still no luck.
The error I get is the following:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/Milevski/.qgis2/python/plugins\ViewshedAnalysis\viewshedanalysis.py", line 194, in run
    z_obs_field, z_target_field, curv, refraction)
  File "C:/Users/Milevski/.qgis2/python/plugins\ViewshedAnalysis\doViewshed.py", line 496, in Viewshed
    mx_dist = numpy.sqrt(temp_x[:,None] + temp_y[None,:])
MemoryError

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.6.0-Brighton Brighton, exported


texts:
Does anyone know a good alternative to the IMapServer2.ExportScaleBar function?

I want to do essentialy the same task, but outside the ArcGIS Server environment. For a given map scale, predefined scale bar style, scale bar size and image resolution, I'd like to generate a scale bar image.

I know I could add a scale bar to the page layout and then export the whole layout, but I do not want to clutter and/or modify the page layout if at all possible.

This is the closest I was able to get, but works only if the resolution used is 96 DPI (screen resolution).

IMapSurround mapSurround = null;
IMap map = GetMap(); // returns the map

var originalScale = map.MapScale;

try
{
    var scaleBar = DefaultSymbols.ScaleBar; // returns a scalebar with some predefined properties
    var scaleBarUid = new UIDClass() {Value = scaleBar.GetType().GUID.ToString("B")};

    mapSurround = map.CreateMapSurround(scaleBarUid, scaleBar);

    var resolution = 300;

    map.MapScale = 1000;

    mapSurround.Refresh();
    ((IScaleBar)mapSurround).UseMapSettings();

    var querySize = (IQuerySize)mapSurround;
    double widthPoints = 0;
    double heightPoints = 0;

    querySize.QuerySize(ref widthPoints, ref heightPoints);

    var widthInches = 3.0; // scale bar is 3 inches wide
    var heightInches = widthInches * (heightPoints / widthPoints); // determine height by respecting the aspect ratio

    var pixelsWidth = (int)(widthInches * resolution);
    var pixelsHeight = (int)(heightInches * resolution);

    using (var bitmap = new Bitmap(pixelsWidth, pixelsHeight, PixelFormat.Format32bppRgb))
    {
        bitmap.SetResolution(resolution, resolution);

        using (var graphics = Graphics.FromImage(bitmap))
        {
            graphics.Clear(Color.White);

            var simpleDisplay = new SimpleDisplayClass();
            simpleDisplay.DisplayTransformation.Resolution = resolution;
            simpleDisplay.DisplayTransformation.Units = esriUnits.esriInches;
            simpleDisplay.DisplayTransformation.Bounds = new EnvelopeClass
                 {
                     XMin = 0,
                     YMin = 0,
                     XMax = widthInches,
                     YMax = heightInches
                 };
            simpleDisplay.DisplayTransformation.VisibleBounds = simpleDisplay.DisplayTransformation.Bounds;

            var tagRect = new tagRECT {left = 0, top = 0, right = pixelsWidth, bottom = pixelsHeight};
            simpleDisplay.DisplayTransformation.set_DeviceFrame(ref tagRect);

            simpleDisplay.StartDrawing(graphics.GetHdc().ToInt32(), 0);

            var drawBounds = new EnvelopeClass
                 {
                     XMin = 0,
                     YMin = 0,
                     XMax = widthInches,
                     YMax = heightInches
                 };

            mapSurround.Draw(simpleDisplay, null, drawBounds);

            simpleDisplay.FinishDrawing();
        }
        bitmap.Save("c:\\image.bmp", ImageFormat.Bmp);
    }
}
finally
{
    map.MapScale = originalScale;

    if (mapSurround != null)
    {
        map.DeleteMapSurround(mapSurround);
    }
}

For resolutions other than 96 DPI this does not work correctly, the scale bar IS drawn to the bitmap, but the measurements it shows are not correct.

Is there a way to do this?


texts: Esri OrthoMaker cut fill volume calculation

texts:
I also posted this to GeoNet as Ortho Maker cannot calculate Cut and Fill.
What I want to do is to calculate a cut fill volume based on my 2 different times DSM which was generated in Ortho Maker as using "Calculate volume in Ortho Maker" following this doc:
https://enterprise.arcgis.com/en/portal/latest/use/cut-and-fill-volume-calculation.htm
Apparently Ortho maker is able to use only their orthoelevation service:
https://elevation3d.arcgis.com/arcgis/rest/services/WorldElevation3D/Terrain3D/ImageServer
Because they let me to use this if I don't have a orthoelevation service in this doc:
https://enterprise.arcgis.com/en/portal/latest/administer/windows/ortho-maker-setup-and-configuration.htm
However, it is silly.
Why would I want to use Esri orthoelevation service to calculate cut fill based on temporal data?
I created Orthomosaic, DTM and DSM for two different times of drone images. I just want to difference from 2 different DSM to calculate cut fill.


texts:
I'm investigating the options available to my employer for the development of a Geological Field Data Collection tool. We would like to develop an Android app that can run disconnected from the internet (we will be operating in some of the remotest places in the world where there are no wireless communications available) and allow teams of field geologists to input their findings into either a surficial or bedrock geological database.

Does anyone have any advice, know of things to stay clear of that could cause problems, know of any software that will help us meet our objective... or any other knowledge to impart that relates to this project.

I will give a summary of what I have come up with so far as well as a bit of background info of the system design that I have envisioned that will work with what we have already.

We have a Windows Mobile ArcPAD application that runs on Getac Mobile PC/GPS and was used previously for the last phase of our project, but being that several years have passed and technology has progressed we want to take advantage of touchscreen capabilities and get away from having to interface with teeny little icons & keyboard with a stylus.

Okay, so here's what I see:

There is a laptop for daily route/traverse planning with a repository of shapefiles to be used as a background/basemap for the system. The team leader opens an application that allows him/her to select an area and draw a "traverse" (linear route for the team to follow) on it and then load this onto a tablet. The leader does this for up to 10 different 2 person teams per day.

The tablets are then distributed to their teams and they head out with the tablets GPS activated and guiding them to the beginning of their traverse.  Once they are there they follow the route laid out for them and identify geological features along the way.  The features are points lines or polygons and have attributes that are available to be chosen from drop down lists, radio buttons and checkboxes. This will standardize the language used to describe the features identified.

At the end of the day the teams return their tablets, and the data management leader then imports and incorporates the days findings into an intermediary geodatabase that holds the data until the crew returns to civilization.  Once they are back the data is loaded into a master Geological spatial database with surficial and bedrock manifolds.

Does anyone know of anything that I should be aware of?

I'm thinking:

Develop a front end program with C# with Xamarin that uses SQLite to store the geological data.

Is there an open source C# shapefile viewer out there that could be incorporated easily?

How hard is it to use the GPS in the tablet to collect data?


texts: Finding empty attribute field with ArcPy search cursor?

texts:
I'm trying to iterate through all fields in my feature classes to find which fields are empty.

I have tried using the SearchCursor by borrowing bits of code from other answers.

This is what I have:

if fieldList:
with arcpy.da.SearchCursor(fc2, fieldList) as cursor:
    for row in cursor:
      if all(row)== NULL:
            print row.name + "is empty"

Where am I going wrong?


texts:
In QGIS 2.18

I have two shapefile polygon layers, the actual case is quite complicated so I'm going to simplify it for the question. One is a layer of large polygons, say counties or states, and the other is a layer of smaller polygons , like individual properties or postcode districts. The polygons in the larger layer have an attribute representing a target, say the desired number of trees on each property. The smaller polygons layer has the actual number of trees each property has. I'd simply like to carry out some process which will highlight those properties/districts which are not meeting the target in their respective counties/states.

There are about 40 such targets I need to compare, so options that involve complicated transformations will not really solve the problem for me. What I need is some sort of expression I can select by which compares the value in one layer to another (I can already do this with the dbvalue() function), but only compare it to the value for the polygon it is contained by.


texts: Query a polyline layer using arcGIS's Javascript API

texts:
I'm attempting to query a polyline layer using the ArcGIS Javascript API I all ready have a lot of experience doing this with Point layers. So I set my Polyline layer up the same way and stopped and started my service again to get every thing going. When I run the query (specifically when I run PolylineQueryTask.execute()) I get an error that just says the operation can't be completed. There is some code below that shows how I'm setting my query layer up. Any ideas on what I can do to fix this issue or things I should try?

PolylineQueryTask = new esri.tasks.QueryTask("http://" + DomainAddress + "/ArcGIS/rest/services/SNND/MapServer/1");

PolylineQuery = new esri.tasks.Query();
PolylineQuery.returnGeometry = true;
PolylineQuery.outFields = ["ObjID"];

PolylineQuery.where = 'where 1 = 1';
PolylineQueryTask.execute(PolylineQuery, ShowPolylineResults, FatlError);

So  as you can see from my code I should be getting all my polylines back. Any thoughts are very appreciated.


texts: Calculate Acres in Territory using ArcGIS Online

texts:
I have a map of counties with each county have assigned a corn acre attribute number. I also have a territory layer that covers multiple counties. I am curious if there is a way to calculate the amount of acres to the territory based on summing the area the territory covers in ArcGis Online.

For Example: County A has 10,000 acres but territory covers 50% of the area of the county then the territory only sums 5,000 acres?


texts: Specify language for mapbox tiles in leaflet

texts:
I want to show the map in my app in English (no localized names in the language most people can't read)

L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token={accessToken}',

I'm using mapbox tiles and it's working for me, but all the signs are shown in local languages. Since my app is in English, I don't want to show Chinese signs to my visitors.

I found the docs about how to change the language in mapbox.js, but I don't want to use their library. Is there a way to build tiles url so it would've been aware about the language preferences?

update: mapbox is not really an important requirement. I'm looking for the tile with labels in English. Can't find any


texts: Network Analyst hierarchy routing issues. (ArcMap)

texts:
My friend and I have created a walking network. We created a network without a hierarchy evaluation and the routing worked fine, but when we added a hierarchy to the equation many routes won't solve. We created a hierarchy class and attributed each of our walk types with a hierarchy of 1, 2, or 3. We thought that the cause of the routes not working was when a route would jump from a hierarchy of 3 to 1, or 1 to 3, skipping 2. But routes do work in some of those situations.

Does anyone know why it breaks, or if there is a way to fix it? (We know that all of our paths are snapped together)


texts:
I have a problem when I add a OpenLayers.Layer.WMS, I get a picture of my map but I can't do zoom on it. I read a lot of questions here but nothing help me. My MapServer map file is:

MAP
 ##CONFIG "MS_ERRORFILE" "/ms4w/tmp/ms_error.txt"
 ##DEBUG 2
 ##CONFIG "CPL_DEBUG" "ON"
 NAME mapa2
 SIZE 1400 900
 IMAGETYPE png
 IMAGECOLOR 196 240 255
 EXTENT -16.1737 27.6847 -13.1254 29.4506
 SHAPEPATH "datos"
 FONTSET "C:/../fonts.txt"
 WEB
      TEMPLATE "/../index.html"
      IMAGEPATH "/../ms_tmp/"
      IMAGEURL "/ms_tmp/"
      METADATA
         "WMS_TITLE" mapa2
         "WMS_SRS" "EPSG:4326 EPSG:4258"
         "WMS_ONLINERESOURCE" "../mapa2.map"
         "WMS_ENABLE_REQUEST" "*"
      END
 END
 PROJECTION
    "init=epsg:4258"
 END
 UNITS DD
 #FOTOS AEREAS
 LAYER
    NAME Foto1
    TILEINDEX indice.shp
    TILEITEM location
    DATA Foto1
    TYPE RASTER
    STATUS on
    OFFSITE 252 252 252
    METADATA
        WMS_NAME "Foto1"
        WMS_FORMAT "IMAGE/PNG"
        WMS_SRS "EPSG:4258"
    END
 END
 LAYER
    NAME    "MUNICIPIO"
    STATUS  ON
    DATA    MUNICIPIO
    LABELITEM "NOM_MUNICI"
    METADATA
        "WMS_TITLE" MUNICIPIO
    END
    PROJECTION
        "init=epsg:4258"
    END
    TYPE POLYGON
    CLASS
        STYLE
            COLOR 254 226 197
            OUTLINECOLOR    255 0 0
        END
        LABEL
            SIZE SMALL
            COLOR 0 0 0
            ALIGN CENTER
            REPEATDISTANCE 9
            POSITION CC
        END
    END
 END
END

And I try to add the layer:

wms = new OpenLayers.Layer.WMS(
  'Municipios',
  "../mapa2.map&service=wms&mode=map&request=getcapabilities",
  {layers:"MUNICIPIO"},{
   isBaseLayer:"true",
   srs: "EPSG:4258",
  format: 'image/png',
  visibility: "true",
  transparent: "true",
  singleTile: true,
  ratio: 1 });

Anybody can help me?


texts: Summing values from different intersecting/overlapping points in QGIS

texts:
I am currently digitizing an old archaeological excavation. I have created the 2m grid where finds were excavated, and have created centroids within each grid. All the corresponding finds has then been joined to the centroids and saved as a new point layer:

The point layer contains multiple points stacked upon each other, because they share the same "local coordinates" ie. 450/560.

What I want to do, is to sum the amount of "something" within between different columns.
Yellow marks "Amount" and where two points share the same coordinates

So I should like to merge all rows with the same X and Y values and have their amount summed.

I should in the end, like to end up with a new polygon layer, where each m2 contains the SUM of all the finds from that specific m2, based on different query's. So it will be possible to graduate colors based on finds intensity


texts: Selectively extracting symbology from layer file in ArcGIS?

texts:
I'd like to classify polygons within a shapefile according to one of my attribute fields, lets call it LITHOLOGY. I have another field in the same shapefile that is associated with a standardised colour code layer file (.lyr), the field is filled with colour codes, lets call this field COLOUR.

I can apply the correct colour symbology to my shapefile by importing the .lay file symbology linked to the COLOUR field. Great. Now, the labels for the symbology are just colour codes. I'd like to apply this colour symbology (which is linked to my COLOUR field to my LITHOLOGY field (which have the right names).

How can I do this?

When importing the symbology from the layer layer file, my options to import "just the symbols or "just the classification" are greyed out.

When trying to produce a map, I either have:

a correctly coloured map with an undesirable series of numbers in my legend, or
I have a legend showing my LITHOLOGY field but my map colours are totally arbitrary.

Unfortunately .lay files can't be opened up and modified and they can't be easily converted to any sensible format either. There's close on 1000 different colour categories, so manual editing of each LITHOLOGY unit is out of the question.

Has anybody run into this issue, because I can't find any solutions on-line?


texts: Esri Query Task service does not display

texts:
I am trying to make my way through some ArcGIS server Javascript API tutorials. My current issue is that when I try to display an Esri.Map, all I see is the "powered by esri" logo in the bottom right. My init is below as is my mapDiv.

function init () {

    map = new esri.Map("mapDiv",{
      basemap: "hybrid",
      center: [-80.94, 33.646],
      zoom: 8
    });

    //Create your query
    var queryTask = new esri.tasks.QueryTask("http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Petroleum/KSFields/MapServer/0");
    //set the onComplete event handler, in this case, when the query is complete, call initLineID,
    //production code would handle handle the error callback as well
    dojo.connect(queryTask, "onComplete", initLineID);

    //build and execute your query
    var query = new esri.tasks.Query();
    query.outFields = ["field_name"];
    query.text = "all";
    query.returnGeometry = true;
    queryTask.execute(query);

}

Any reason why this would not be working?


texts: PostGIS: Single side ST_Buffer of linestring produces strange output

texts:
In the example below I have a linestring PostGIS table (in black and yellow) and a single sided ST_Buffer output polygon (in red). The buffer distance is based on an attribute value stored in the linestring table.
As you can see one of the buffers is incorrect, and cuts a section of itself out for a certain amount of it. Why has this happened? Can it be fixed? Seems to be an issue where the linestring section becomes concave.

Try this for yourself using a similar query:
SELECT ST_BUFFER('LineString (638327 331059, 638337 331048, 638337 331037, 638343 331023, 638358 331014, 638380 330991, 638382 330978, 638391 330966, 638410 330959, 638420 330953)',50,'side=right')

In another example I'm seeing this happen (the black line is 1 feature, as is the buffer polygon):

And try this one with this:
SELECT ST_BUFFER('LineString (600296 172961, 600304 172962, 600307 172965, 600307 172969, 600307 172971, 600309 172974, 600308 172976, 600307 172979, 600306 172981, 600307 172982, 600311 172980, 600313 172979, 600317 172979, 600321 172982, 600325 172982, 600327 172979, 600328 172976, 600331 172974, 600338 172974)',10,'side=right')

Increase the buffer distances for stranger results. EPSG: 27700 in both cases.


texts: PyQGIS: looping over shapefile's fields to rename them using substring

texts:
I'm trying to rename multiple field names (about 150) of a shapefile in QGIS. Their names are different but they all end by _xf20 and I would like to remove this part for each of them.
I've looked up to other topics that could help and I've found this one: Rename fields of shapefile using PyQGIS 3 with the code:

layer = iface.activeLayer()

# Open editing session
layer.startEditing()

# Rename field
for field in layer.fields():
    if field.name() == 'oldName':
        idx = layer.fields().indexFromName(field.name())
        layer.renameAttribute(idx, 'newName')

# Close editing session and save changes
layer.commitChanges()

and this one: Delete all fields with the string "RAWI" using QGIS Python
But even if it is a start, I'm a beginner with Python and I don't find a way to adapt them for my case.
I would rather not edit every single field manually.


texts: Which points are decisive in the "Oriented Minimum Bounding Box" algorithm in QGIS?

texts:
I used the ombb algorithm in QGIS and there were no problems with it. I am just wondering which points determine / influence the angle, how the orientation of the box is determined and when a box is inclined by more than 90 ° or less than 90 °. The documentation doesn't help me.
In ArcGIS there is the function "CalculatePolygonMainAngle" and I think that the angle of the ombb can vary a lot from the ArcGIS function. I don't have ArcGIS and unfortunately I can't try it out.


texts: How to reverse vertex ordering in polygons?

texts:
I am working with an Oracle spatial connection in QGIS through which I am accessing and editing a vector layer. I find that when I edit a geometry, the vertex ordering is reversed (to clockwise in the outer ring) which then gets flagged up as an error in oracle.

I am using the combine function (http://www.qgis.org/api/classQgsGeometry.html) among others in my own editing tools, but I notice that this happens with the native merge polygons tool too.

I think this reo-ordering to clockwise is the default in QGIS based on this link: https://hub.qgis.org/issues/6283 and was wondering if there is a way to reverse it to counter clockwise? (preferably with python)


texts: Uploading older Google Earth Historical Imagery into QGIS

texts: There are numerous ways to import Google Earth imagery into the QGIS (e.g. OpenLayer plug-in or XYZ tiles), however the recent-most images are downloaded. Is there any way to import historical imagery? My project's related to tracking land use dynamics and it would be great if I could switch between imagery from different years.


texts: Rasterise Polygons with consistent resolution and size in QGIS

texts:
Using QGIS I want to rasterise a series of polygon layers so the output is the same size and resolution.  I am using Rasterise, however am not getting the result I want.

I have a raster that is the size and resolution that I want - "Raster".

I have used Rasterise with the options
Input file - Polygon layer I want to rasterise "MMEIL"
Attribute field "Sus_Lvl"
Output file for rasterized vectors "Raster"
Keep existing raster size and resolution
Load into canvas when finished.

The problem is this updates the raster in the areas where I have the polygon, but the rest keeps the values of the previous raster.  I need to overwrite the entire raster.  My Sus_Lvl has values 1,2,3.  I want a raster that has as output values in the range 0 (no data), 1, 2,3.

If I simply set the width and height I get what I want in terms of output - but I cannot keep a consistent resolution between all the different layers.

A solution would be to update the existing raster so that all values are 0 - but I cannot figure out how to do this.


texts:
The following picture shows my routing database all set for querying loaded to QGIS. Red line is the road I want to select route from, yellow points are the vertices generated after running pgr_createTopology() and the violet points are bus stops. As per the tutorial here, everything works fine using qgis db manager.

I want to show route from one bus stops to another using pgRouting. Since there are more bus stops than the automatic generated vertices, one vertices is selected for multiple bus stops.Is this is because of something wrong in my road shapefile or anything else? is there anything wrong with the procedure i am following ?


texts: Using value from same row in Field Calculator of QGIS

texts:
I am using QGIS.
I have an attributetable with several columns.
I want to change the values of column1 if they are NULL with the values of the same row in column5.
case WHEN column1 is NULL THEN
END

How do I say "use the value from the same row"?


texts: Adding Combobox to Python Toolbox tool dialog?

texts:
I am trying to create a Python Toolbox tool which is supposed to add two fields State and StateAbb to selected Shapefile from selected combobox.

How I can add a Combobox which contains a list of US states to the Tool?

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "states"
        self.alias = "states"

        # List of tool classes associated with this toolbox
        self.tools = [states]

class states(object):
    def __init__(self):
        self.label = "states Tool"
        self.description = "Add and Insert Value from Combo"
        self.canRunInBackground = False

    def getParameterInfo(self):

        param0 = arcpy.Parameter(
            displayName="1- State Street Layer",
            name="in_features",
            datatype="GPFeatureLayer",
            parameterType="Optional",
            direction="Input"
        )

        params = [param0]

        return params

    def execute(self, params, messages):
        """The source code of the tool."""
        arcpy.AddField_management(params[0],"area","TEXT","#","#","#","#","NULLABLE","NON_REQUIRED","#")
        arcpy.CalculateField_management(params[0],"area","!shape.area@squaremeters!","PYTHON_10.3","#")

        return


texts:
Is there a way that i can upload my CartoDB map in Openlayers as a WMS layer or Vector layer. Here is the my CartoDB map link.
http://cdb.io/1nxdrpf


texts: PyQgis: During editing, add new feature attributes

texts:
Workflow: Start edit mode -> New feature is created -> On the console (pyqgis), new attributes are set -> Exit edit mode.

Problem: The attributes are not saved. Only after completion of the editing, I can edit the attributes.

The feature is created and during the editing, I try the following:

newFeature = layer.selectedFeatures()[0]
 layer.dataProvider().changeAttributeValues({ newFeature.id() : {0 : 'foo'} })
 layer.updateFields()

The Layer already has an attribute table and some features.
After editing, the above code works. During editing, no attributes can be set.

Before the feature appears in the map, it opens a window. In this window, attributes can be entered. This also works. But I would like to change the attributes later, after the feature has been displayed in qgis.


texts: Problem with Euclidean distance in ArcGIS

texts:
In ArcGIS 10.7.1, when I run the Euclidean distance it completes but the result is blank and doesn't show anything in data view and the distance in table of content is 0.
When running the Euclidean tool the environment setting were:
Extend: as default
Coordinate system: same as the layer
How can I solve the problem?


texts:
I have a dozen raster (3 bands) files in GeoTIFF format (32bit), which I have downloaded from SentinelHub (sentinel2 satellite imagery).
I try to compress them before I publish them as an Imager Service through ArcGIS Server. I follow the guidelines from this link, which makes use of the "Translate (Convert Format)" tool.
The output raster is very small (500kb) and it is all black. I thought I need to calculate statistics for it, but this didnt fix the issue.
Any ideas what could be the issue?


texts: Getting script file name as object in Python?

texts:
I'm quite new to Python and programming in general.

I have created a few scripts using Python and as part of my framework/best practices for my scripts, I'd like to display when a script begins and ends in the IDLE shell.  For now, I have just been adding a print function with the name of the script.  So for example for a script called "UpdateMaps.py", I would have:

import arcpy, os

print "BEGIN: 'UpdateMaps'"

<some script that updates maps>

print   "COMPLETED: 'UpdateMaps'"

But I would prefer to have the script name be printed automatically, since I often cut and paste pieces of script to another one, change the name of the scripts, etc.

So I would prefer something like this (I'm making up the 'str(scriptname)' object as an idea of what I'm looking for):

import arcpy, os

print "BEGIN: " + str(scriptname)

<some script that updates maps>

print   "COMPLETED: " + str(scriptname)

Is there a way to do that?


texts: Calculate variable importance for classifications in Google Earth Engine

texts:
I am doing a simple Random Forest classification in Google Earth Engine. I  try to classify upland, water and wetlands, using 25 variables derived from one year of Sentinel-1 images. The classifier works well and the accuracy is okay, but now I want to calculate a variable importance for my classifier.

There seems to be no function in Earth Engine for calculating variable importance. The only workaround I could think of is exporting the stacked image to R, and do the classification there. However, this might be a complicated solution, due to the size of the study area and the high number of variables. Is there a simple work-around to calculate variable importance in Earth Engine?


texts:
According to Opening MPK file without ArcGIS Desktop?, .mpk files are just .zip compressed and I should be able to unzip them as normal. I've tried from the terminal on macOS (10.11.6):

unzip somefile.mpk

for multiple different .mpk files I always get this error:

Archive:  somefile.mpk
 End-of-central-directory signature not found.  Either this file is not
 a zipfile, or it constitutes one disk of a multi-part archive.  In the
 latter case the central directory and zipfile comment will be found on
 the last disk(s) of this archive.
unzip:  cannot find zipfile directory in one of somefile.mpk or
        somefile.mpk.zip, and cannot find somefile.mpk.ZIP, period.

I get the same error if I change the file extension from .mpk to .zip. These files are known to work fine on Windows, when uncompressed with ArcGIS tools. They can also be unzipped with 7-zip on Windows.

Any idea what's going on here?


texts:
I'm doing a project that evaluates the practicality of using remote sensing imagery in agricultural field management.

I noticed there are some scratches/strips on PlanetScope imagery for one of the paddocks. I'm not sure what causes this and more importantly, is it fixable and how? If not fixable, how can I avoid downloading bad quality imagery?

Product: 20180108_051245_0c79_3B


texts:
I have over 17,000 individual polygons. I want to merge two adjacent polygons so that I finish with approximately 8,500 polygons. The polygons need to be merged based on location, merge two adjacent polygons. Then merge the next two adjacent polygons and so on. Is there a way to do this through ArcGIS or Python?

This is a sample of the python code I am trying to use

#Import arcpy
import arcpy,os,datetime,re
import math
from arcpy import env

#The filepath location
mxdPathA3P="C:\\Users\\heff\\Desktop\\ArcGIS_data\\HSCU_Mapping\\HSCU_A3_Portrait_Template.mxd"

#Set the workspace environment
arcpy.env.workspace = 'C:\\Users\\heffernanv\\Desktop\\ArcGIS_data\\ARCPY_WORK'

inpgs=arcpy.mapping.ListLayers(mapdoc,u'POLYGON_LAYER',mainDF)[0]

# tempMerged = temporary merged polygon shapefile (better to have this in-memory)
where1 = '"' + 'Shape_Area' + '" < 125000000.0000'
shapeField = arcpy.Describe(inpgs).ShapeFieldName

#Read polygon geometry into dictionary; key = GridID, value = geometry
geometryDictionary = {}
polyrows = arcpy.SearchCursor(inpgs)
for prow in polyrows:
    geometryDictionary[prow.GridID] = prow.shapeField
del prow
del polyrows

#Loop through catchments that are to be dissolved
polyrows = arcpy.SearchCursor(inpgs, where1)
for prow in polyrows:

    #Get geometries of two polygons
    # grd1 = GridID on the current polygon
    # grd2 = GridID on the polygon downstream
    geometry1 = prow.shapeField
    geometry2 = geometryDictionary[grd2]

    #Merge datasets
    arcpy.Merge_management([geometry1, geometry2], tempMerged)

    #Dissolve
    arcpy.Dissolve_management(tempMerged, outpgs)


texts: Algorithm to find geohashes which satisfy query condition

texts:
I have an application where I need to query geohashes by specifying the bounding rectangle. Since geohashes map 2D to 1D is there a method by which I may find the prefixes of all geohashes which satisfy the query condition (bounding rectangle). I tried to google but I did not find any algorithm.

It will be great if there is some library which does that for me efficiently in C/c++. Even if there is no library then too I can can work with an algorithmic solution.


texts: How to import a background image for overlay with GeoPandas

texts:
I recently downloaded the shapefiles for San Francisco from Mapzen. Now I am able to display the streets, but I would like to have a real map as a background image, such that it is easier to make the link between the real roads from San Francisco and the ones shown by GeoPandas.

Here is the code to show the streets with GeoPandas:

First download the shapefile from Mapzen: https://s3.amazonaws.com/metro-extracts.mapzen.com/san-francisco_california.imposm-shapefiles.zip
Open a Python notebook and add those lines:

 #Useful starting lines
 %matplotlib inline
 import numpy as np
 import matplotlib.pyplot as plt
 %load_ext autoreload
 %autoreload 2

Load the street map:

shapefile_dir = 'the path to your directory where you store the shapefile'

shapefile_name ='san-francisco_california_osm_roads_gen0.shp'
shapefile_roads = os.path.join(shapefile_dir, shapefile_name)

Import as a geopandas dataframe:

import os
import geopandas as gpd

df = gpd.GeoDataFrame.from_file(shapefile_roads)

Simply plot the street network:

df.plot();

The output should be:

Now I would like to have a real map (Color map) as a background.
How can I do it?


texts: Would PostGIS offer an advantage over MySQL for a produce farm application?

texts:
I have a web app that stores the locations of farms in West Michigan. You can search for a product (e.g. "broccoli") and it will show you all the farms that grow that product.

Right now I'm using MySQL and using trigonometry to calculate the difference between the user's location and the location of each farm. It's not a bad way to go but it did take some doing.

Another thing I want to do soon is to map out growing seasons for different products for different regions. (For example, I want to show that avocados grow at a certain time of year in California but never Ohio.)

I realize this is an open-ended and possibly naive question, but might it be worth it for me to make the switch to PostgreSQL/PostGIS to take advantage of its spatial capabilities?


texts: Unit testing with postgis 1.5 and django

texts: I'm trying to get some unit tests working in Django/PostGis. When the test database is formed, it doesn't have the postgis and postgis_topology extensions installed. Is there a way to issue commands to the database before the creation of tables happens? Or could you somehow mess with the postgres template to inherit from? Or something else that wouldn't require me to either hack the django code or not destroy the database? Thanks!


texts:
My professor gave us a basic Python script that demonstrated how Shapely works with a shapefile, however it's not working for me. Here is the code:

from shapely.geometry import mapping, shape
from fiona import collection

with collection("hawaii_vector_data.shp","r") as input:
  for polygon in input:
    c = polygon['properties']['class']
    if c == 0 or c == None:
      continue
    s = shape(polygon['geometry'])
    print s

and when I run it I get this Traceback error:

Traceback (most recent call last):
  File "test.py", line 6, in <module>
    c = polygon['properties']['class']
KeyError: 'class'

Since I am new to using Shapely I have no idea what the error might be.  Any advice would be great.


texts: Spatial Database Design (PostGIS) - Point data that doesn't always have spatial data

texts:
I am creating a spatial database using PostGIS.

I am storing the attributes and location of some species records that can be represented as points. The records are recorded from either:

a) static survey position (where previously I have just represented the species point as the static position) or
b) a transect walked around the site and a position is marked when a record is recorded.

The database stores the static positions in a a separate table, so for a) above I don't actually have to store a point against each of these records as I can do a join on the species table and the static location table. If I did store the location, if the static location changed, I would have to also make sure I moved the location of all of the species records associated with that point.

The transect results are stored with spatial data as they could be recorded anywhere around the walked route and there is no other way to know where they were recorded.

I would all of this data in the same table as it seems silly to have to maintain two tables for records about the same species (would have to make changes to the table structures in two tables).

I can see three ways to deal with this:

Store the species data in two separate tables, one spatial and one not.
I am guessing you can't have NULL spatial values in a geometry column?
Give all the records spatial data, and move the points from the species table when the static table locations change? Can you have a relationship / procedure in the database that automatically does this for you?

Which of these approaches is the most sensible? Or have I missed something. Any thoughts would be much appreciated.


texts: Draw middle ordinate perpendicular to chord, out to curve, with arcpy

texts:
I've been assigned to find the radius of road curves. My boss found this article describing a pretty straightforward method of doing so:

find the start/end points of the curve;
create a chord;
create the middle ordinate (from chord midpoint, perpendicular line that goes up to the curve)
do a bit of math with the chord and ordinate lengths to find the radius.

Manually, this works very well. The downside is that I'm working with thousands of miles of road, so we want to build an automated tool to be able to quickly analyze our entire existing network, and also be able to check roads digitized in future.

I have adapted this script to find transition points where the curve begins/ends -- comparing line segments to see if the road direction is changing (e.g. is it still turning right, or has it started curving left instead). I therefore have chords for all the curves, but now I'm struggling with how to make middle ordinates. It needs to:

start at the midpoint of the chord
be perpendicular to the chord (I'm not positive how to do this yet, but found a script that will probably work and/or be adaptable)
be long enough to reach the road, but no longer (???)

So, long question short: how can I find the (approximate) distance between my chord midpoint and the road?

(Another StackExchange question describes working with the Near_analysis function which looks like it could work, but I currently have all of those chord midpoints just stored as a python list, not a feature class. I'd like to avoid generating a bunch of new features if possible.)


texts: gdal_merge.py cannot import gdal in Windows 7 Anaconda GDAL Installation

texts:
I have an installation of GDAL through Anaconda that is mostly working, but it throws an ImportError whenever I run gdal_merge.py:

Traceback (most recent call last):
  File "C:\Users\xxx\AppData\Local\Continuum\Anaconda2\Scripts\gdal_merge.py", line 37, in <module>
    from osgeo import gdal
ImportError: No module named osgeo

I can import GDAL without error from my own scripts and commands, even within the above folder. Most of what I have read about similar issues is that they are usually related to the PATH environment variable, but as far as I can tell mine is pointing to all the right places:

C:\Users\xxx\AppData\Local\Continuum\Anaconda2;C:\Users\xxx\AppData\Local\Continuum\Anaconda2\Scripts;C:\Users\xxx\AppData\Local\Continuum\Anaconda2\Library\bin;

I have 64-bit Anaconda for Python 2.7 installed on a Windows 7 system. My version of GDAL is 2.1.0, which I installed with the conda install utility. I set and tested my GDAL_DATA environment variable, and believe it to be correct:

C:\Users\xxx\AppData\Local\Continuum\Anaconda2\Library\share\gdal

I have replaced my user name with xxx throughout.


texts: I use two different Landsat datasets in different computers. But the region settings I work is the same. Is it possible to export the region settings in the form of a file in GRASS?


texts:
I am using the WGS 84/UTM coordinate sysytem.

When I go to > Project Properties / General / Canvas Units to Degrees, Minutes, Seconds my coordinates are wrong.

I am using QGIS 2.0.1


texts: Saving shapefile with QgsVectorFileWriter

texts:
I've tried to create a vector file (shapefile) with QgsVectorFileWriter. This is my code:
path = "/home/camilocorredor/Documentos/ANT/4_Abril/File.shp"
options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "ESRI Shapefile"

QgsVectorFileWriter.writeAsVectorFormatV2(layer_test, path, QgsCoordinateTransformContext(), options)

The vector file never is created on the folder.
Any suggesting?


texts:
I try to add MS SQL Server Tables to QGIS on Ubuntu 18.04.3 LTS

I installed msodbcsql17 and msodbcsql17

I configured the .ini files

odbcinst.ini:

[SQL Server]

Description=Microsoft ODBC Driver 17 for SQL Server

Driver=/opt/microsoft/msodbcsql17/lib64/libmsodbcsql-17.4.so.1.1

UsageCount=1

and odbc.ini :

[kplus]

Driver =SQL Server
Server = ip-adress, 1433
Database=kplus

I can connect with isql command isql -v kplus pdw usr and query tables.

When I add MS SQL Connection I receive the error :


texts: FME assign specific node id with OSM writer

texts:
Is there a way to reassign the node's id value to something I generate?

I did try exposing the uid attribute in the write transformer however only produced a new attribute and did not affect the id value.

//The question has been changed. To clarify the question at hand, How can I/Can I assign a node's id with a specific value when fme starts to write an osm file?

Thanks.


texts: If the Earth was flat, what unit of measurement would distances between geographic points be represented as?

texts: A whole lot of abstraction- given two longitude/latitude points given from a GPS, I'm using the Pythagorean distance formula to find the distance between them. What unit of length would the distance be? Kilometers? Miles?


texts: r.to.vect is not responding in QGIS GRASS

texts:
NLCD data entered for half of Illinois, all of Michigan and Wisconsin, plus additional along borders.
Used r.neighbors to simplify, then r.to.vect with "smooth corners" feature selected.
Got one warning

Concurrent mapset locking is not supported on Windows.

Got to "Registering primitives..." and has been there nearly three days.
Task Manager indicates 25% to 30% CPU utilization, approximately 500 MB or more Memory use with Power usage High to Very High. Running on i3-2379M CPU @2.4 GHz with 6.00 GB Memory, 64-bit Windows 10.
Should I let it run or "End Task"?


texts: Composer map needs to reset extents each time

texts: I have several maps that are regularly updated and am now having a problem that has only recently begun - each time I reopen a composer to re-PDF my maps, I have to reset the map extents.  All of my maps have been created the same way and I've never had any issue before, but now there's several of them where the map frames will just be white and not display any data until I reset to the map canvas extents.


texts: How does GeoServer/RasterSymbolizer ramp interpolate?

texts:
I want to encode data using a colormap using GeoServer/RasterSymbolizer, by mapping data using the full 24 bits in RGB. I am thinking of the following scheme:

<ColorMap type="ramp">
    <ColorMapEntry color="#7FFE00" quantity="-512"/>
    <ColorMapEntry color="#7FFFFF" quantity="-256"/>
    <ColorMapEntry color="#800000" quantity="0"/>
    <ColorMapEntry color="#8001FF" quantity="256"/>
    <ColorMapEntry color="#800200" quantity="512"/>
    <ColorMapEntry color="#8003FF" quantity="768"/>
</ColorMap>

However, I am not completely sure this will work, as I could not find documentation on the RasterSymolizer. Between two color values, only the B channel (last two symbols) should change.

A value of 257 using this colormap, would be interpolated between #8001FF and #800200. I hope this would encode to #8001FE, but theoretically it could also be #8002FE.

Alternatively, I could double the entries in the colormap to make it more explicit, like so:

<ColorMapEntry color="#7FFEFF" quantity="-256"/>
<ColorMapEntry color="#7FFF00" quantity="-256"/>
<ColorMapEntry color="#7FFFFF" quantity="0"/>
<ColorMapEntry color="#800000" quantity="0"/>
<ColorMapEntry color="#8000FF" quantity="256"/>
<ColorMapEntry color="#800100" quantity="256"/>
<ColorMapEntry color="#8001FF" quantity="512"/>
<ColorMapEntry color="#800200" quantity="512"/>

Is it allowed to have two entries for the same quantity?

A third option is like below:

<ColorMapEntry color="#7FFEFF" quantity="-257"/>
<ColorMapEntry color="#7FFF00" quantity="-256"/>
<ColorMapEntry color="#7FFFFF" quantity="-1"/>
<ColorMapEntry color="#800000" quantity="0"/>
<ColorMapEntry color="#8000FF" quantity="255"/>
<ColorMapEntry color="#800100" quantity="256"/>
<ColorMapEntry color="#8001FF" quantity="511"/>
<ColorMapEntry color="#800200" quantity="512"/>

How will a value of 255.5 will interpolated in this case? #800180 or #800080? Following the encoding scheme it should be #8000FF. Note that the input data is stored as a geotiff with float32's.


texts: Leaflet Markercluster and Timeline in sync

texts:
So I have been trying to figure how to use the markercluster and timeline plug in in sync to show a temporal animation of points. As some of the point data overlaps a cluster is needed to be able to select all the data. Unfortunately I have been having issues with using the pointToLayer function to create a marker and add to my marker cluster.here is a sample of what I tried with no luck:

    var timeline = L.timeline(data, {
      formatDate: function(date){
        return moment(date).format("YYYY");},
     pointToLayer: function(feature, latlng){
      var marker = L.marker(latlng);
      marker.bindPopup("test");
      markers.addLayer(marker);
      return markers;
     }


texts:
I understand there’s a way or a program that can assign polygons (precincts) to a feature point.(Vote Centers)

For insistence the law states we need to assign a maximum of 10 precincts (or polygons) to a Vote Center. We have 607 precincts and 69 Vote Centers.

I utilized the "join" but it in some instances it assigned more than 10 polygons(precincts) per point(Vote Center)
Where do I set parameters (minimum or maximums)?


texts:
So I've tried using the slope tool for a DEM using an appropriate Z-factor and the output is smooth, but it seems to be over-estimating slope angle.  I would prefer to reproject the DEM (to NAD 83 UTM Zone 9), and run the same analysis, but the output has line artifacts:

I've tried calculating statistics and using bilinear interpolation, but to no avail.


texts:
I try to create some advanced expression using a QGIS plugin for education rules.

I have a simple Gui :

My Gui has 7 line edits for text inputs.

I try to using this Expression:

layer = QgsMapLayerRegistry.instance().mapLayersByName("facility")[0]
request = QgsFeatureRequest().setFilterExpression('"name" =\'{0}\' and "name1" =\'{0}\' and "name2" =\'{0}\' and "name3" =\'{0}\' and "name4" =\'{0}\' and "name5" =\'{0}\' and "name6" =\'{0}\' '.format(lineedit1,lineedit2,lineedit3,lineedit4,lineedit5,lineedit6,lineedit7))
for i in layer.getFeatures(request):
        print I

But if some line edit (it is not necessary to have text in all line edits) is empty then my expression doesn't return some result.

If I try something like this :

if lineedit1 !='' and lineedit2=='' and lineedit3==''and lineedit4==''and lineedit5==''and lineedit6==''and lineedit7=='' :
    layer = QgsMapLayerRegistry.instance().mapLayersByName("facility")[0]
    request = QgsFeatureRequest().setFilterExpression('"name" =\'{0}\''.format(lineedit1)
    for i in layer.getFeatures(request):
            print I

or :

if lineedit1 !='' and lineedit2!='' and lineedit3==''and lineedit4==''and lineedit5==''and lineedit6==''and lineedit7=='' :
    layer = QgsMapLayerRegistry.instance().mapLayersByName("facility")[0]
    request = QgsFeatureRequest().setFilterExpression('"name" =\'{0}\' and "name" =\'{0}\''.format(lineedit1,lineedit2)
    for i in layer.getFeatures(request):
            print I

.....................

then the expression works fine. But I have 7 inputs and many combinations if I continue this example.

Any idea how this expression could work?


texts: Are there any simple map data sources like world.dat

texts:
I recently came across this blog post www.gnuplotting.org/plotting-the-world/ that used a file that holds a basic map of the world. Is there a name for this kind of data. Also is there a repository of similar but higher resolution data.

Direct link to data:
www.gnuplotting.org/data/world.dat


texts: Bonn 3.2.2 to buggy need to install older version of QGIS

texts:
I am having a lot of issues with Bonn 3.2.2 QGIS and would like to uninstall this version and install an older version that is more stable.

What seems to be a stable version of version 3?
What is the best way to do so on Linux?


texts: PyQGIS: feature added signal throws an error message

texts:
For my PyQGIS plugin 2.2, I try to:

def logFeatureAdded(self,layerId,addfeatures):
    message = str(layerId) + " has features added: "
    for feature in addedFeatures:
        message += str( feature.id() ) + ", "
    QgsMessageLog.logMessage(message)
    QApplication.beep()

    layer.featureAdded.connect(self.logFeatureAdded)

When I try to add a feature, this following error message appears:

  TypeError : logFeaturesAdded() takes exactly 3 arguments (2 given).

How can I use the featureAdded signal properly?


texts:
I want to create combobox with all available scales, that could replace standard zoom slider control. As far as I know, there is no scale attribute in OL and scales must be calculated from resolutions. In api docs I found that resolutions can be retrieved from:

view - if fixed resolutions were provided during initialization. In my case they were not and the result is undefined

var resolutions = this.map.getView().getResolutions();

tileGrid of a layer - I do not have any base layer and for each layer in my map code below does not return resolutions

resolutions = layers[0].getSource().getTileGrid().getResolutions();

So now, how can get all available resolutions? I believe I can calculate them from projection but I can't find any proper method to do that.

EDIT:

Immediately after posting the question I found a workaround. I can create TileGrid from projecton and get the resolutions from there:

var projection = ol.proj.get(this.getSrs());
        var tileGrid = ol.tilegrid.createXYZ({
            extent: projection.getExtent(),
            tileSize: 512
        });

        console.log('projection resolutions', tileGrid.getResolutions());

but it does not seem to be good solution, so maybe somone can point me to the right direction.


texts: How to filter the legend according to the current Atlas page?

texts:
I filter the legend by map in the map composer. Some symbols must be used for two legend entries, as you can see below.

When I start the atlas tool, the legend displays - if the symbol contained by the map - both legend entries. Is it possible to filter the legend by the atlas and not by the map.


texts: How to publish many MXD files in ArcGIS Server 10.1 automatically?

texts:
How to publish many mxd files in ArcGIS Server 10.1 automatically?

I have about 60 mxd files and I want to publish them automatically.
Is there any way to do this?


texts:
I want to get the centroids of my multipolygons using Python and GeoPandas. But when viewing the polygons and the centroids points on QGIS, the centroids are not inside the polygons.
Here is the code I used:
import geopandas as gpd
import pandas as pd
import sys

# Read the data
polygons = gpd.GeoDataFrame.from_file('mypolygons.shp')
temp=polygons.centroid
polygons['x'] = temp.x
polygons['y'] = temp.y
centroid=polygons[['x','y']]

gdf = gpd.GeoDataFrame(centroid, geometry=gpd.points_from_xy(centroid.x, centroid.y))
gdf.to_file("centroid/centroid.shp")

What is wrong?


texts: Importing Feature Class to Oracle SDE Geodatabase Simplifies Features

texts:
I recently imported a feature class from a File Geodatabase to an Oracle SDE Geodatabase. The features of the two feature classes now look drastically different despite not having any edits made.

The green polygon is the original feature as it looked in the file geodatabase. The purple transparent polygon on top is the feature after being imported into the Oracle SDE Geodatabase. It's as if the act of importing the feature class simplified the features.

Any ideas on why that would happen? I am new to Oracle and ArcSDE, maybe I set something up incorrectly?

Here is another view of the same area, with different symbology:

The storage format is ST_GEOMETRY, and there was no conversion performed, as far as I know, I simply imported the feature class into the Oracle SDE geodatabase from the file geodatabase.


texts:
As far as I understand, the m value "is a measured distance along a route" (similar answers here), thus similar to the mileage of a car driving over a highway. The m value defines the point on a line as the distance from a set starting point. So the m-value, even if not explicitely set, is somehow inherent in the very geometry of the line feature, defined as the distance from the first vertex of the geometry. Thus, every point of the line (not only vertices) has a well defined m value - comparable to an x and y coordinate, even if these are explicitely defined only for the vertices, not the points on the segments between vertices.
However, with a look at the QGIS documentation, I am confused. The image to explain the Filter vertices by M value function there shows a line with 5 vertices, with m-values in this order: 0, 11, 5, 12, 8. This would correspond to the mileage going up and down, so something is wrong.
How is this possible if the route is measured from the start point of the line? Or did I misunterstand the concept and the m value relates only to the neighboring vertex, not to the start point of the geometry?


texts: Input & near features fields into the point distance table

texts:
wondering in ArcGIS 10.3x how can I carry-over the fields in the "Input Features" and "Near Features" layer to the point distance table? The point distance table only lists the Input_FID, Near_FID, Distance.

In above, the point distance table only has FID and distance but the "subs" and "lld" fields are missing...


texts: QGIS PostgreSQL connection working but very slow on mobile connection

texts:
I'm working with a QGIS project while commuting with train, it works rather alright. The problem is that rendering a PostgreSQL table on my mobile connection takes very long time compared when I'm on my home wifi. (I have a PostgreSQL server at home that I store data in)

I usually have a rather fast mobile connection (50Mbit/s down and 5Mbit/s upp) but I guess that I have a rather high loss. Is there any setting or something I can do to improve the rendering speed?


texts: Performing specific table joins in QGIS via python?

texts:
I am trying to update an existing field on a shapefile using data from another shapefile. Basically, each attribute table contains an ID field that I want to use to update a field called "address".

So far, I have tried to use the processing alg "qgis:joinattributestable" and QgsVectorJoinInfo(). Both of these solutions however seem to create temporary join fields that are of no use to me.

Is this something that is possible with QGIS?


texts: Fails to run JoinAttributes in stand-alone PyQGIS script

texts:
I'm writing a stand-alone script using no GUI. I cannot use the function Join Attributes by Location. I think that the problem might be with proper passing arguments to the function. There's a piece of the code:

from qgis.core import *
from qgis.utils import *
import os,sys

QgsApplication.setPrefixPath("'C:/OSGEO4~1/apps/qgis'", True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append('C:/OSGEO4~1/apps/qgis/python/plugins')
print sys.path

from processing.algs.qgis.JoinAttributes import *
alg = JoinAttributes()

alg.setParameterValue('INPUT_LAYER',QgsVectorLayer("C:\Points.shp",'pkt','ogr')) #all of those layers exist
alg.setParameterValue('INPUT_LAYER_2',QgsVectorLayer("C:\Polygons.shp",'polygons','ogr'))
alg.setOutputValue('OUTPUT_LAYER', QgsVectorLayer("C:\Ewid.shp",'CadastralLayer','ogr'))

from processing.gui.SilentProgress import SilentProgress
progress = SilentProgress()
alg.processAlgorithm(progress)

The error is as follows:

Original exception was:
Traceback (most recent call last):
File "mainScript.py", line 33, in <module> aCD.addCadastralData(layers,parcelLayerPath)
File "Z:\Soft\AutomatyczneRaportowanie\AddingCadastralData\AddingCadastralData_main.py", line 36, in addCadastralData alg.processAlgorithm(progress)
File "C:/OSGEO4~1/apps/qgis/python/plugins\processing\algs\qgis\JoinAttributes.py", line 76, in processAlgorithm layer.crs())
File "C:/OSGEO4~1/apps/qgis/python/plugins\processing\core\outputs.py", line 308, in getVectorWriter crs, options)
File "C:/OSGEO4~1/apps/qgis/python/plugins\processing\tools\vector.py", line 550, in __init__ if self.destination.startswith(self.MEMORY_LAYER_PREFIX):
AttributeError: 'NoneType' object has no attribute 'startswith'

Could you help me to sort it out?


texts: ArcGIS Select Analysis using Unique Values from Field

texts: I am trying to create separate shapefile using Select analysis with ArcPy by doing an SQL query to create them based on a unique value from the field. I want to create a list of the unique values in the table from a specific field since there are some duplicates of the same value in the field. I want to do a search of the table and find all the unique values and then append them to the list and then do the Select analysis to create individual shapefiles and eventually individual KML files. Thanks in advance.


texts: Defining one numeric variable by defining another string variable (via parameter selection)

texts:
I'm creating an ArcGIS tool that will allow users to figure out how much nitrogen that a field needs based on interpolated points (Kriging). I will have the following parameters:

#numeric variables
yieldGoal = arcpy.GetParameterAsText(0)
organicMatter = arcpy.GetParameterAsText(1)

#text variable
previousCrop = arcpy.GetParameterAsText(2)

To figure out how much nitrogen is already in the soil, I have to do some Kriging:

n_krig = Kriging(soilPoints, field, kModelOrdinary, cellSize, kRadius, outRaster)

The next step is to calculate the amount of nitrogen that the soil needs. The equation takes into account previous crops, because some can add nitrogen to the soil, and others do not.

I want users to be able to choose the previous crop that was grown with the previousCrop = arcpy.GetParameterAsText(2) parameter, and based on that selection, I want to create another numeric variable pca that will serve as a nitrogen amount for that previous crop selection. The 'pca' variable can then be plugged into the final equation:

fertilizer = ((1.6 * yieldGoal) - (organicMatter * 20) - (n_krig) - (pca))

How can I define pca?

I tried to write an ifstatement to define the variable, but it doesn't work:

#Original parameter:
previousCrop = arcpy.GetParameterAsText(0)

#If statement
if previousCrop == str("Corn"):
    pca == int(0)
elif previousCrop == str("Wheat"):
    pca == int(0)
elif previousCrop == str("Milo"):
    pca == int(0)
elif previousCrop == str("Sunflowers"):
    pca == int(0)
elif previousCrop == str("Soybeans"):
    pca == int(40)
elif previousCrop == str("Alfalfa: >5 plants/square foot"):
    pca == int(120)
elif previousCrop == str("Alfalfa: 2-5 plants/square foot"):
    pca == int(80)
elif previousCrop == str("Alfalf: 2-3 plants/square foot"):
    pca == int(40)
elif previousCrop == str("Alfalfa: <2 plants/square foot"):
    pca == int(0)
else:
    pca == int(0)


texts: Calculating geometry for a clipped layer

texts:
Working with ArcGIS Pro. I'm trying to calculate stream length for areas outside of state land. I've tried this two separate ways and was successful in getting the layer that I'm looking for, but for some reason, I go to calculate geometry and it spits out the length for the full stream not the modified/short stream (I know this since I can compare IDs and their respective lengths). This seems like it should be simple... Am I missing something? How do I calculate the new stream length?
Here are the two ways that I have created my working layers:

I pulled in my stream layer and used the Erase tool to remove the streams that overlapped the state land layer.

I used the Clip tool with the stream layer to the state land layer. I figured if I could calculate the stream length from the streams within the state land, then I could substract that from the total stream length and get stream length outside of the state land.