texts: Spatial Join in PostGIS according to Line Attribute

texts:
I am in the process of building a data warehouse in PostGIS, that will store the length of the road network of a country (Lines) per Administrative Area (i.e. Province or State, referred to as AA hereafter and is a Polygon). This is my first time of building a large scale Relational data warehouse so I am not really familiar with this kind of operations.

The "id" field of the spatial join resulting table is used to relate to the unique "id" field of the AA polygon layer. The length of the road segments is already calculated and is derived from the "len" attribute of the Road Network table.
So far the operation yields the desired results by performed the following query:

SELECT
  m.id,
  m.order08,
  sum(r.len) as length,
FROM
  public.and_lines_2012_03 AS r,
  public.and_a8_2010_12 AS m
WHERE
  ST_Intersects(m.geom,r.geom)
GROUP BY
  m.id, m.order08
ORDER BY
  id;

And the resulting table looks something like this:

gid |   id   | length
 1  | 000002 | 74118
 2  | 000012 | 128131
 3  | 000013 | 174296
 4  | 000014 | 82240
 ...|  ...   |  ...

The Road Network is divided in 3 different line types that are attributed in the Road Network table by the name "line_type" and can have values 0,1,2
My intention is to create a table that will have a column for the length of each line type for each Administrative Area. That is the table would ideally look something like this:

gid |   id   | length_type_00 | length_type_01 | length_type_02
 1  | 000002 |     30028      |     11564      |     32526
 2  | 000012 |     64270      |     15602      |     48259
 3  | 000013 |     104492     |     51981      |     17823
 4  | 000014 |     38153      |     9524       |     34563
 ...|   ...  |      ...       |      ...       |      ...

I have experimented with different approaches in SQL statements, but so far have not managed to grasp the approach of this kind of queries.
The most relative resource I have found on the Internet is this one from OpenGeo:

Update 2012/05/23

@Nicklas

I have managed to produce what seems to be very accurate results with the way you have suggested. The query used is the following:

SELECT
  b.id AS id,
  b.order08 AS order08,
  round(sum(b.line_type_00)) AS line_type_00,
  round(sum(b.line_type_01)) AS line_type_01,
  round(sum(b.line_type_02)) AS line_type_02
FROM
(SELECT
  m.id,
  m.order08,
  (line_type=0)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_00,
  (line_type=1)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_01,
  (line_type=2)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_02
FROM
  public.and_lines_2012_03 AS r,
  public.and_a8_2010_12 AS m
WHERE
  ST_Intersects(m.geom,r.geom)
GROUP BY
  m.id, m.order08,line_type
ORDER BY
  m.id
)
AS
  b
GROUP BY
  b.id,b.order08
;

The query is quite slow though (4846 tuples take around 41 seconds), since from what I understand it re-projects on the fly to produce the the results in meters. Have I understood/am I doing something wrong?


texts:
I need to test a condition within the sql in postgis.
My sql:

SELECT
    (zoneamento3_32722.zt3),
    logradouros32722."NMIDEABR",
    logradouros32722."OBJECTID",
    "public".logradouros32722."NRIMPINI",
    "public".logradouros32722."NRIMPFIN",
    "public".logradouros32722."NRPARINI",
    "public".logradouros32722."NRPARFIN"
FROM
    "public".logradouros32722,
    "public".zoneamento3_32722
WHERE
    st_within (
        logradouros32722.geom,
        zoneamento3_32722.geom
    )
AND logradouros32722."NMIDELOG" LIKE 'IPIRANGA%'

/ / if odd goes here
AND (
    logradouros32722."NRIMPINI" <= 60
    AND logradouros32722."NRIMPFIN" >= 60

/ / if pair goes here

AND
 (
        logradouros32722."NRPARINI" <= 60
        AND logradouros32722."NRPARFIN" >= 60
    )
)


texts: Earth Engine invalid band id when exporting multiband image

texts:
I am trying to export a multiband image to my assets using this code:

var aoi= ee.Geometry.Rectangle([-1.2, 41.9, 2.8, 43.2]);
// filter 1 Jan to 30 Jan
var c = ee.ImageCollection("MODIS/006/MOD10A1")
      .filterBounds(aoi)
      .select('NDSI_Snow_Cover')
      .filter(ee.Filter.calendarRange(1, 30, 'day_of_year'));

// export collection to multiband image
var stack = c.toBands().clip(aoi);
Export.image.toAsset(stack, 'stack', 'stack', null, null, null, 500);

The Task details window returns the error message: Error: Invalid band id '2001_01_01_NDSI_Snow_Cover' specified.

I have tried many things like renaming the bands, changing the data type, etc. without success. Another post reported a similar issue but I am not using the colon character here..


texts: Unit of earth radius while calculating distance based KNN weight matrix using PySAL

texts:
I have geographical polygon data in WGS84 (lat/long) [EPSG: 4326] geographic coordinate system. I want to calculate distance-based KNN weight matrix using pysal which I can calculate using following code.
from libpysal import  weights
w_knn = weights.distance.KNN()

But there is a radius parameter which corresponds to the earth's radius. pysal uses this to remove effect of earth curvature while calculating distance between geographic units.
But in the official documentation there is nothing about the unit of this radius. Can anyone suggests something about unit of this radius?
I also want to know that if our dataset is projected in a projected coordinate system (say UTM) then should I ignore this radius parameter while calculating weight matrix.


texts: Calculating mean DNI (Direct Normal Irradiance) for bird species in QGIS

texts:
I am new to GIS so excuse my lack of proper terminology.

What I would like to do is to calculate a mean DNI (direct normal irradiance) for multiple bird species. I do find a DNI-raster online (https://solargis.com/maps-and-gis-data/download/world) and have a multipolygon (?) containing species-presence polygons (available from birdlife). I've added a picture of the layers to ensure that my terminology is correct.

I would now like to create an attribute that is basically the mean irradiance a bird species would receive given the DNI that its distribution range receives. What I envision is that this would mean that each "pixel" the species occurs in should be multiplied by its corresponding DNI value. Next all these "pixels" should be summed and then divided by the total amount of occurence-"pixels" for that species. I am just not sure how to do this.

I hope the question is clear, and that someone could help me with this.


texts: Automate publishing ArcMap mxd to ArcGIS Online using ArcGIS REST API or Python

texts:
I have a mxd file that contains several feature classes and they get updated every day. I am with ArcGIS 10.1 at the moment. It would be potentially upgraded to 10.5 in the near future anyway.

I have successfully published the mxd file to ArcGIS Online as a Feature Layer using ArcMap -> Share As -> Service -> Hosted Services on ArcMap.

As the feature classes get updated every day, I would like to achieve an automated approach to allow the mxd file to be published/updated to ArcGIS Online on a daily basis. Such that the Feature Layer can be updated there accordingly.

Is this something that can be achieved using ArcGIS Python or REST API?


texts:
I have some Python3 code used to import shapefiles using fiona. This is the main set of lines:
I created some code (Python3) with these lines:
c = fiona.open(dataDir + "/" + sourceFile, encoding=shapefileEncoding)
for i in range (0, len(c)):
                        # Get next record
                        record = next(c)
                        .......
c.closed

The code works perfectly but the output includes this warning:
FionaDeprecationWarning: Collection.__next__() is buggy and will be removed in Fiona 2.0. Switch to `next(iter(collection))`.  record = next(c)

I probably don't understand what iterable iter refers to. I replaced next(c) in my code with next(iter(c)) and my code no longer works: I get the same record repeatedly instead of iterating through all records.
What should I replace the next(c) in my code with?


texts: How to define bbox value in TileCache tool?

texts:
I have created a WMS layer on GeoServer and i would like to tile it using TileCache tool, now. There are some definitions about bbox parameter in the web page of the tool but i can't understand.

In GeoServer, there is bbox like depicted below; (btw: projection is 3857)

And the configuration of TileCache is;

[LayerName]
type=WMS
url=http://address:8080/geoserver/wms?layers=LayerName&TRANSPARENT=TRUE&tiled=false
extension=png
srs=EPSG:3857
size=256,256
levels=22
bbox= ?
maxResolution= ?
extent_type=loose

I don't know how can i define bbox and maxResolution values. There is one more problem; tiling bbox value! If i would like to run command below, what should be the bbox value again?

> tilecache_seed.py LayerName 1 3 -b "(?) 3873263.280,4494756.762,4119659.272,4809662.117" -p 1


texts: Mosaic multiple 7gb TIFF into one big TIFF

texts:
I have data of 10 sets each with 361 channels in format of ENVI.
Each set contain:

7gb format-less file (BSQ #UPDATE)
2.4gb .ENP file
9kb .HDR file

I try to mosaic all of them in one big file.
What I tried:

Used ENVI to mosaic all of them but it takes forever...
Used different parameters to speed up the process but no help...
Tried rasterio, worked on small test set (used converted to tif files via ENVI).
But it doesn't work with the full set, I think GDAL can't handle such big files.
ArcGIS Pro, the result end up messed up channels (the 361)


texts: QGIS 3.10 attribute table change to default or modify when export the layer to different format

texts:
When I export my newly digitalized layer I noticed, that my data attribute changes.
In the QGIS I set the data attribute table like this:

but each single export to the different format change it, as it appear as a new separate layer in QGIS.

For instance, when I exported my file as a .kml layer, then my data attribute table looks like this:

My digitalized layer with set attribute layer is already active (1), all other are exports to the different formats, where this data attribute table have changed (2)

For example in Google Earth it looks like this:

Is it something, which I should do between these steps?


texts: How to enter series of points along certain angle at regular distances in ArcMap?

texts: I would like to enter a series of points along a certain angle in ArcMap.  I have a UTM coordinate as a starting point, and from that I'd like to enter 10 points at a 120 degree bearing every 50 meters. What is the best way to do this?


texts: Is it possible to consume ArcGIS Hosted Table into local or ArcSDE geodatabase?

texts:
I would like to publish a table to our ArcGIS Online (AGOL) organization account.

Is it possible to consume hosted features in AGOL to a local db or ArcSDE database connection?

By consume I mean import a hosted ArcGIS Online table into a geodatabase with ArcPy.

By local db, I mean file geodatabase


texts: How to enable smooth zoom in Openlayers in default library set?

texts:
I've compared two Openlayers samples:

Online dev example - clientzoom
release dev example - clientzoom

First one zooms in smoothly, second one just flashes and shows tiles from next zoom level.

The only difference between those two versions is the library - release dev example (second) uses default multi-file version, while the first one - uses some compiled version (not sure from which tag).

Question is: how to enable smooth zoom in Openlayers 2.12? Or how to enable smooth zooming in any Openlayers version ;)

[Edit]

The fluent zoom feature is included by default in 2.13 library:

release dev example for 2.13 - clientzoom


texts:
I'm looking to automate a procedure that will change a single cell value in a raster.  I have a known x, y and z, and a raster.  I need to change the cell value (closest to my x and y value) to the specified z value.

I've been doing this using Saga (change cell values [interactive] module) but I need a more automated method for use in a batch file or python script.
MI thought there may be a module for updating a raster cell based on point (vector) values?

Note, I do not want to reclassify grid values, rather I wish to specify one.

I have access to QGiS, Saga and ArcGIS basic (without any spatial or 3D packages).


texts: I have 33 layers, 1 with points and 32 with lines. I want to extract/highlight the points where more than one line is intersecting. I have tried with buffers but I'm not getting anywhere.


texts: MapBasic: skip popup menu dialog to get access to standard "Modify table structure" dialog

texts:
I'm trying to develop a small program using MapBasic for MapInfo.
I want to get access to the standard "Modify table structure" (or "Image registration" for raster tables) dialog through "Layer Control" menu using right click on the layer. It may be very helpful especially when you've got a lot of open layers and tables in MapInfo.

Standard way is a little bit challenging:

You need to Remember the name of your layer
then you go to "Table" -> "Maintenance" -> "Table structure..." -> Find your table in the popup menu ->
Begin to modify it.

The same thing with the raster registration: you always must to find your table in the popup menu before you can modify it.

I want to skip this popup menu.  You just right click on the layer, choose "Modify table/Raster Registration" option and go straight to the "Modify table structure" dialog. The issue is that I cannot figure out how to get access to these standard MapINfo dialogs. I have found them in "mires
.dll" using "resourse hacker", but I still cannot skip the popup menu.

Is it possible?


texts: Get selection geometry from select interaction

texts:
I have a select interaction that is drawn in the map as a polygon while the map feauters was not loaded because of scale limit on the layer.

When I zoom in so that the features load, I want to apply a selection style on the newly loaded features that is inside the selection polygon.

I have a custom loader on the data layer, so I can get the selection interaction during the load. The problem is that I can't find any api to get the polygon geometry from the selection interaction. I do have the interaction object available.

How can I solve this?


texts:
I'm trying to load around 6000 shapefiles (all the postcodes for a city) into GeoDa but it crashes everytime ("OpenGeoDa.exe has stopped working"). Their troubleshooting page details this problem but I've not been able to resolve the issue, as stated here: http://geodacenter.asu.edu/node/404#load

I've tried to code a new id field in Arcmap but that did not work (was I right in recoding the new ID variable as a short integer?).

I then tried in Arcmap to see if the field names were the problem. I cleaned up most of them up, making them shorter and removing underscores. However, I have been unable to change 5 of the fields (could they be read only?) and the names are over 30 characters long.

Could this be the problem?

Please see the screen grab below for an overview of my data and the length of one of the field names.

Screen grab number 2. Note: I reduced the number of fields to 5 and only selected a few polygons to see if size was the issue. The GeoDa error message on the right is displayed (still the same message).


texts: Why postgis st_overlaps reports Qgis' "avoid intersections" generated polygon as overlapping with others?

texts:
I am implementing a postgis trigger on polygon layer to check and reject newly created polygon if it overlaps others.

The problem shows when I digitize new polygons adjacent to others in Quantum GIS with option "avoid intersections" in Settings -> Snapping Options:

Sometimes Postgis allow me to save and some time it doesn't with a
message saying "overlapping polygon detected" (programmed into
trigger).
Sometimes Postgis doesn't allow saving the first time but
when I save again (same unmodified geometry), it's done without any
notice.

My trigger function:

CREATE OR REPLACE FUNCTION check_overlapping() RETURNS TRIGGER AS
$BODY$
DECLARE
BEGIN
IF TG_OP = 'INSERT'
THEN
  IF
    (SELECT COUNT(*)
     FROM
       (SELECT id
        FROM polygon_layer AS t
        WHERE st_overlaps(NEW.geom, t.geom)) AS foo) > 0
  THEN
    RAISE EXCEPTION 'OPERATION ABORTED: overlapping polygons!';
    RETURN NULL;
  ELSE
    RETURN NEW;
  END IF;
ELSIF TG_OP = 'UPDATE'
THEN
  IF
    (SELECT COUNT(*)
     FROM
       (SELECT id
        FROM polygon_layer AS t
        WHERE st_overlaps(NEW.geom, t.geom)
        AND (t.id <> OLD.id)) AS foo) > 0
  THEN
    RAISE EXCEPTION 'OPERATION ABORTED: overlapping polygons!';
    RETURN NULL;
  ELSE RETURN NEW;
  END IF;
END IF;
END;
$BODY$
LANGUAGE plpgsql;

CREATE TRIGGER check_overlapping_trg
BEFORE INSERT OR UPDATE ON polygon_layer
FOR EACH ROW EXECUTE PROCEDURE check_overlapping();

My software version: Qgis 1.8.0, Postgis 2.0

Added June 24th, 2013:

Create table query:

CREATE TABLE polygon_layer
(
  id serial primary key,
  geom geometry(MultiPolygon,3857)
);
CREATE INDEX polygon_layer_idx
  ON polygon_layer
  USING gist
  (geom);

Test data:

insert into polygon_layer (geom) select st_GeomFromEwkt('SRID=3857;MULTIPOLYGON(((11805973.8551485 1136129.8216546,11757217.3210366 1107981.71948687,11790391.87002 1102117.53153526,11823733.9672305 1116526.6790735,11801114.95656 1119207.45070852,11805973.8551485 1136129.8216546)))');
insert into polygon_layer (geom) select st_GeomFromEwkt('SRID=3857;MULTIPOLYGON(((11819880.3580052 1063581.43928182,11791397.1593831 1067602.59673435,11799483.2424957 1095075.79477947,11809995.0126011 1081844.19604541,11843504.6580388 1127249.76561359,11817272.092427 1129138.51033764,11832111.37859 1146852.90819468,11881203.0091563 1142329.10606058,11819880.3580052 1063581.43928182)))');
insert into polygon_layer (geom) select st_GeomFromEwkt('SRID=3857;MULTIPOLYGON(((11809995.0126011 1081844.19604541,11792995.1839403 1103242.5817722,11823733.9672305 1116526.6790735,11801114.95656 1119207.45070852,11804235.9486147 1130077.11269212,11843504.6580388 1127249.76561359,11809995.0126011 1081844.19604541)))');

Trigger behavior:

When run 3 insert queries, the first 2 will go smoothly but the 3rd will raise exception as defined in trigger function.

If I drop the trigger, insert all 3 polygons and check to find overlapping polygons with following query, I can not find anything:

select *
from polygon_layer as a
     inner join
     polygon_layer as b
     on st_overlaps(a.geom, b.geom)
where a.id < b.id


texts: I am fairly new to ArcGIS and just wanted to know if anyone knew of any particularly helpful tutorials online or books for making story maps?


texts: Making spatial clusters by direct manipulation of raster pixel values using R

texts:
I've processed a raster layer with pixel values ranging between 1 & 8 (see here). A cursory glance reveals that, there's a good tendency of clustering in different parts of the heterogeneous raster band, occasionally intruded by unwanted unlike pixel values. My goal is to make spatial clusters with the following characteristics-

Every spatial cluster will represent only one pixel value i.e. all the minor intruding pixel values will be converted to the prominent pixel value surrounding it;
There can be more than one cluster with the same pixel value, but has to be spatially well-separated;
The clusters can be of any size &/or shape.

Here you can see a few clustering example zones drawn by hand to show tentatively how I want them to be classified. I have searched the 'raster' package of R, but failed to discover anything relevant. A painstaking solution to such problems is digitizing, which, in many ways, is unprofessional, especially when someone has to deal with similar situations over and over again.


texts: How do I report area statistics for land cover classes at equal sections along a river?

texts:
I have a river outlined:

Using QGIS, the Semi-automatic Classification plugin and Landsat 8 images I classified the habitats (blue=water, yellow=sand/rock...,green=vegetation/algae) inside the river polygon across multiple years:

I then created 50 sections along the length of the river dividing it equally:

How do I derive habitat area statistics in meters (using WGS84 CRS for the layers) for each of these equal sections?

I want to relate habitats to river bird abundances, comparing each sample's (river section's) habitat to its birds!


texts: Error message: "Image.rename: The number of names must match the number of bands. Have 4, found 1."

texts:
I want to estimate the VCI (Vegetation Condition Index) but I get the error message:

  Image.rename: The number of names must match the number of bands.  Have 4, found 1.

Any ideas? I think, the problem is at the final part of the script: "VCI CREATION"

Here is the Code Editor script: https://code.earthengine.google.com/008283891289f6d518058b2d9f2175c4


texts: Does raster to polygon conversion impact area calculation using calculate geometry?

texts:
I am a beginner GIS user.
I am trying to calculate the area of clipped buffers which I have intersected using the intersect tool by a habitat raster which I am working with. The edge of the intersected buffers just hangs over the edge of the raster (see image below). I used the clip (data management) tool for this and then I converted the output clipped raster file into a polygon using the raster to polygon tool. However, now when I try and calculate the area of the clipped intersected buffers it seems to be very different from what I was expecting. Is this due to the raster to polygon conversion? Any suggestions about how I can get around this?


texts: Portal for ArcGIS - Group Administrator only

texts:
Would it be possible to have Portal for ArcGIS local administrators being able to manage their group only?

Basically, admin of group A should not be able to view or edit group B. I want to assign administrator level 2 role only for the group he/she owns.  The Portal administrator will create private groups and assign admins for each groups to manage their own content and membership.


texts: Using Specific MapServer with ArcGIS for Server?

texts:
I'm trying to use this MapServer, but I'm having a very difficult time with the attribute table of a specific layer. I need the coordinates of each hydrant in franklin county in kentucky.

When connecting to the server through ArcGIS and opening the hydrant layer all I see are the hydrants. I need the coordinates of each hydrant. I know the data is there, because when I click on one of the points using the identify tool the information that should be in an attribute shows.

Map Server: http://agis21.frankfort.ky.gov/ArcGIS/rest/services/Hydrants/MapServer/0

I've tried using the Metadata toolset with no luck.


texts: I work with line vector data (road network).The lines (roads) includes several vertexes. I would like to have for each line only nodes - at the end of line, or at the intersection of n lines. How can I delete all the vertices, to have clean lines with nodes only, but keep the topology? When I delete vertex, the line topology is changed. My QGIS version is 3.8.2.


texts: How to get from a 7-folder-deep directory structure to a 3-folder structure for TMS?

texts:
so I got a Tiled Map Service up and running using TileCache to store all my tiles and OpenLayers to actually show them in a browser using an Apache server. But the demo map tiles I found were stored in a 7 folder deep directory, e.g.

13/000/002/123/000/004/936.jpg

would be zoom level 13, column 2123, row 4936. But I've been given tiles formatted just 3 folders deep, e.g.

13/2123/4936.jpg.

I'm guessing the first, 7-segment is the TMS standard since it worked correctly and everything displayed correctly, but I was wondering if there was any easy fix to display the 3-segment directory correctly, or am I just going to have to change it to fit the 7-folder format? And is there any easy way to do that? Because that'd be a lot of files to have to change. Thanks.


texts:
Using the below I output a KML file from a Large SpatialPolygonsDataFrame

polygonWGS <- spTransform(Polygon, CRS("+proj=longlat +ellps=WGS84
+datum=WGS84"))
writeOGR(polygonWGS, dsn="polygonWGS1.kml", layer="Area", driver="KML",
overwrite = TRUE)

However, the resulting attributes table is blank in this KML. Specicially I'm interested in having the polygons named as the 'Area' name.

Is this possible?


texts: Clear of vector is not happening in OpenLayers

texts:
I have a screen where I have a button 'A' on click of which a popup opens. In that popup, I have a map which consist a vector layer (Polygon). I have another button 'Remove Shape' in that popup. Now when I click that Remove Shape button, I am able to remove the vector layer from the map by using
this.map.getLayers()"array_".getSource().clear();

Although, the layer (Shape Polygon) is getting removed, on closing the popup and reopening it, I can't remove the shape again.
Can any one help me with this and let me know if I am missing anything anywhere ?


texts: VectorLayer Feature setStyle makes features disappear in OpenLayers 4

texts:
After I have already displayed my map I want to highlight certain features.
For starters I am trying to highlight all features like this:

public highlightMapFeatures() {
    this.map.getLayers().forEach(function (lyr) {
      if (lyr.type === 'VECTOR') {

        lyr.getSource().forEachFeature(function(feature) {
          let style = new ol.style.Style({
            stroke: new ol.style.Stroke({ color: '#ff00ff',
              width: 2 }),
          });
          feature.setStyle(style);
      }
    });

But what happens is, that all my features disappear.
Do I have to redraw them? I should be so simple.


texts: How do I merge shapefiles with picture attachment tables

texts:
I am trying to clean up some GIS files from some much older projects where data collection was done differently. currently we create a shapefile with domains to do inspections, in the past we just used individual shapefiles for certain locations. Each of these past style shp files would have 7 files, I would like to merge these into 1 file (sounds easy enough right) but each of these files has an associated attach table with pictures which need to merge as well (All in the exported GDB)... is there anyway to do this?
I can use either ArcGIS Pro, Arcmap 10.7.1, or AGOL to do this.


texts: Scaling factor for atmospheric correction in GEOMATICA

texts: Wondering if there's any one who's familiar with using ATCOR for atmospheric correction in GEOMATICA. What is the best scale factor? I used the default setting and then divided my image by 1000. This seems to work for one of my images, but when I tried the same value for another image from another area, spectra from corrected image is no where near my in-situ reference data.


texts: Error running TopologicalBreachBurn in Whitebox GAT

texts:
I am trying to execute the TopologicalBreachBurn script in Whitebox GAT. This is my first time using this software so am generally unfamiliar with how it operates. I imported a DEM as an ASCII .txt to convert it to the required .dep format needed for the script. Then I selected a streams vector .shp. When I hit run, it gives me the error: "an error occurred during operations. See log for details."

When I open the log there is nothing there, so I am not able to ascertain what the error actually is. Anyone run into this issue or know what may be happening?


texts: For loop to Create lines of raster mosaics in ArcMap 10.5

texts:
For loop to Create lines of raster mosaics
With the tool clip raster, I am using a for loop code (code included) which works well.

>>> import arcpy
>>> in_raster_folder = "D:\KKR\ct_igeoref"
>>> in_polygon_folder = "D:\KKR\ct_plgs"
>>> out_clip_folder = "D:\KKR\ct_iclipped"
>>> arcpy.env.workspace = in_raster_folder
>>> in_rasters_list = arcpy.ListRasters("*", "ALL")
>>> for raster in in_rasters_list:
...     raster_name = raster [:-4]
...     clip_poly = in_polygon_folder+"\\"+raster_name+"_iplg.shp"
...     out_raster = out_clip_folder+"\\"+raster_name+"iclp.img"
...     arcpy.Clip_management(raster,"" , out_raster, clip_poly, "256", "NONE", "NO_MAINTAIN_EXTENT")

Now using the tool MosaicToNewRaster  I need to create 6 lines of mosaics using each time sets o 6 rasters (basic code included).

arcpy.MosaicToNewRaster_management(input_rasters="D:/KKR/ct_iclipped/018_cticlp.img;D:/KKR/ct_iclipped/019_cticlp.img;D:/KKR/ct_iclipped/020_cticlp.img;D:/KKR/ct_iclipped/021_cticlp.img;D:/KKR/ct_iclipped/022_cticlp.img;D:/KKR/ct_iclipped/023_cticlp.img",
    output_location="D:/KKR/ct_moslines",
    raster_dataset_name_with_extension="ct_line1.img",
    coordinate_system_for_the_raster="", pixel_type="8_BIT_UNSIGNED",
    cellsize="", number_of_bands="3", mosaic_method="LAST",
    mosaic_colormap_mode="FIRST")

My mosaic lines are:

Line1 with 018_cticlp.img + 019_cticlp.img + 020_cticlp.img + 021_cticlp.img + 022_cticlp.img + 023_cticlp.img
Line2 with 026_cticlp.img + 027_cticlp.img + 028_cticlp.img + 029_cticlp.img + 030_cticlp.img + 031_cticlp.img
Line3 with 034_cticlp.img + 035_cticlp.img + 036_cticlp.img + 037_cticlp.img + 038_cticlp.img + 039_cticlp.img
Line4 with 042_cticlp.img + 043_cticlp.img + 044_cticlp.img + 045_cticlp.img + 046_cticlp.img + 047_cticlp.img
Line5 with 050_cticlp.img + 051_cticlp.img + 052_cticlp.img + 053_cticlp.img + 054_cticlp.img + 055_cticlp.img
Line6 with 058_cticlp.img + 059_cticlp.img + 060_cticlp.img + 061_cticlp.img + 062_cticlp.img + 063_cticlp.img


texts: Cannot use Advanced Digitizing Toolbar (greyed out)

texts:
I am new to QGIS. I have one picture with the attribute table I created and the greyed out toolbar. Yes, I have it in "Toggle editing" but doesn't seem to make a difference....even after I close out of the attribute table. Am I missing something?

The shapefile (trap_location) is for points, and agri_landscape is for polygons.


texts: Zonal smoothing on a specified area of a raster?

texts:
Is there any way to smooth a specify an area on a raster, such as a DEM,  and smooth it?
For example, I'd like to smooth out the artifact highlighted in the red box.
What is the best course of action to remove this artifact in ArcGIS? I'm open to using another software if it's free.


texts: Setting up localhost in proxy service for ArcGIS Javascript API 3.x

texts:
I have a proxy service (DotNet) that works great in my production environment but I can't use it in my development environment. Obviously I really need it to work in my development environment because I don't want to deploy something that doesn't work. I get 404 errors for the resource. This is how my current code sits, but I have tried a couple different things which I will explain in a minute:

<ProxyConfig allowedReferers="*"
         mustMatch="false"
         logFile="proxyLog.txt"
         logLevel="Warning">
    <serverUrls>
      <serverUrl url="https://www.mydomain.com/myserver/rest/services/HabitatManagement/HabitatManagement/GPServer"
                 username="xxx"
                 password="xxx"
                 matchAll="true"/>
      <serverUrl url="https://www.mydomain.com/myserver/rest/services/HabitatMonitoring/HabitatData/MapServer"
                 username="xxx"
                 password="xxx"
                 matchAll="true" />
      <serverUrl url="https://www.mydomain.com/myserver/rest/services/HabitatMonitoring/HabitatClassification/GPServer"
                 username="xxx"
                 password="xxx"
                 matchAll="true" />
    </serverUrls>
</ProxyConfig>

I've tried setting the allowed referrer to my specific local host like this:

<ProxyConfig allowedReferers="https://www.mydomain.com/*,"http:/localhost:14906/*"

I have also tried this:

<ProxyConfig allowedReferers="https://www.mydomain.com/*,"http:/localhost:14906/appName/*"

But that doesn't work either. I'm not really sure where to go from here.

This is an example of one of the referring URLs:

  http://localhost:14906/appName/proxy/proxy.ashx?https://www.mydomain.com/myserver/rest/services/HabitatMonitoring/HabitatData/MapServer/2?f=json


texts: Getting compilation error on esriLicenseProductCode.esriLicenseProductCodeArcView

texts:
I am getting compilation error on esriLicenseProductCode.esriLicenseProductCodeArcView

I have installed Arcgis 10.1 and authorized ArcEditor. I am not able to authorized ArcView.
I have not installed Runtime Engine.

Am I missing anything any suggestions.


texts: Merging layers to intersect and count (tally) multiple fields in QGIS?

texts:
My question is similar to Performing multiple layer intersection in QGIS?

However, there are a few differences. I have run a test analysis on a dummy data set and cannot replicate the result with my real larger data set.

I have a list of environmental features that are potentially at risk (streams, fish, amphibians, red listed species, ungulate winter range, etc) that I have merged into a single omnibus layer. I have another polygon study area that consists of the threat (roads and storm water outlets - high salt application area). My goal is to study the intersection between the at-risk-layer and the threat-layer and to obtain a count of the number of at-risk-features that intersect with the threat. In the dummy example everything works just fine (see figure above). The next image is from my real data set.

This looks more complicated. The brownish layer is the merged polygon of streams, wellsites, amphibians, fish in streams, agricultural lands, ungulate winter range, etc... The light grey layer is the study area (roads and stormwater outfalls). The RskCnt field is a sum of all the fields with string data attached. I created multiple fields by running the following:

CASE
   WHEN "ALR" LIKE '%' THEN 1
   ELSE 0
END

Then I added all the fields to get the count:

WellID_n + GrndH2O_n + redlist_n + bluelist_n + wetland_n + riparian_n + uwr_n + Amphib_n + Hydro_n + Fish_n + ALR_n

However, I know that the counts are wrong. There are places where I should have more than a count of 3, because I can look at the individual layers that were not merged and determine that the count is wrong. The merging is not capturing all the fields. The fields are there, but not all the information is getting passed along.

What could be going wrong here?

The difference is that I didn't merge layers in the dummy case, so that's a clue. There must be something going wrong with the merge...but I can't seem to figure this out.


texts: Sentinel-2 active fire detection service

texts:
I'm looking for a Sentinel-2 analogue of the FIRMS service (Fire Information for Resource Management System), which provides the location of active hotspots. The low spatial resolution of MODIS and VIRS data is not enough for me.
The ability to download the data is required.
A quick search on the Internet returned no results.


texts: Transposing coordinates from PDF map onto satellite image

texts:
I have a map in a PDF file with boundaries that I want to plot onto a satellite image.

I have made a screenshot, saved it as tiff and imported it as a Rraster object.
Using click I have outlined the boundaries which I can now plot onto the raster object.

My goal is to convert the data.frame with x and y coordinates to a spatial object (sf or sp) and plot it onto a satellite image.

The problem is the projection: How do I convert my relative positions to a projection that matches that of my satellite image?


texts:
How to give the access rights to access "localhost / lm" since I have the error 403 "forbidden you do not have ....."

Can the problem come from there??


texts: Simple land classification (urban, suburban, rural) using polygons in MapInfo?

texts:
I have a large number of polygons which I need to classify as either rural, suburban or urban. These polygons represent wireless coverage for a cellular network. What I would like to achieve is a breakdown for each of these polygons as to how much rural, suburban and urban coverage they provide e.g polygon1 Rural: 74%, Urban: 0%, Suburban: 26%

Currently my polygons are in MapInfo tab files.


texts: GRASS v.transform tool producing no results QGIS

texts:
I'm working on a simple polygon feature representing the Western Cape province of South Africa. My goal is to enlarge the polygon so it can stand out from the others.

To do so I've used the v.transform tool, but the output layer is empty. The picture below shows how I've computed the tool, I might be doing something wrong.

As I understand it:

X shift: I've put the x coordinate of the centroid of the polygon.

Y shift: same for the y coordinate

X scale: the factor by which I'd like to enlarge my polygon (I've used 130% for the example, according to another post: How to resize shapefile features proportionally in QGIS?)

GRASS GIS 7 region: my canvas extent.


texts: Using QGIS 2.0, I have created a color legend for raster values in the print composer. By default each color swatch appears to the left of its label in the legend but I want the color swatchs to appear to the right of their labels. I'd also like to be able to specify how the items are aligned with the legend's frame (e.g. right/left). I've looked in the item properties for the raster legend in the print composer, but haven't found any controls to effect this behavior.


texts: How to insert polyline on ArcGIS JavaScript 3.x?

texts:
I am using JavaScript 3.x library of ArcGIS for plotting the polyline but nothing to show on the map.

 function simtest() {
                 var p1 = esri.geometry.geographicToWebMercator(new esri.geometry.Point(76.331549,10.013566, new SpatialReference({
                                wkid: 4326
                            })));
var p2 = esri.geometry.geographicToWebMercator(new esri.geometry.Point(76.331463,10.013566, new SpatialReference({
                                wkid: 4326
                            })));
line = new esri.geometry.Polyline();
line.addPath([p1, p2]);
polygraphlayer.add(new esri.Graphic(
  line,

  new esri.symbol.SimpleLineSymbol(
                    esri.symbol.SimpleLineSymbol.STYLE_SOLID,
                    new dojo.Color([124,252,0, 0.1]), 2)
));
var pt = new Point(76.331549,10.013566, new SpatialReference({
                                wkid: 4326
                            }));

 map.centerAndZoom(pt,13);

                }


texts: Change the actual coordinates in a shapefile (QGIS)

texts:
I'm trying to map data for a county in Texas. The data were provided to me in a shapefile, and I've been working with it in QGIS 2.14. My problem is this: The coordinates in the shapefile I was given are wrong. The latitudes should all be in the 31-degrees range, but instead they're around 21 degrees N, which - when exported into Google Earth - puts the county squarely in Mexico.

Is there a way in QGIS that I can simply add ten degrees to all the inaccurate latitude coordinates? Everything else about the shapefile seems correct - it's just the degrees latitude are ten degrees lower than they should be.


texts:
I'm trying to plot GOES-16 data with world boundaries. Data comes from GOES-16 files such as this.
Data comes with "GOES-R ABI fixed grid projection" reference system, which to my better knowledge is describe as this PROJ4 string: '+proj=geos +h=35786023.0 +a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs'
My first step is to create a raster layer:
library(raster)
rrqpe <- raster("OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc", varname='RRQPE')

I get some warnings about not being able to create a valid CRS, perhaps due software version.

Warning messages:
1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
long_name=GOES-R ABI fixed grid projection
perspective_point_height=35786023
semi_minor_axis=6356752.31414
sweep_angle_axis=x
2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
grid_mapping_name; false_easting; false_northing; scale_factor_at_projection_origin;
scale_factor_at_central_meridian; standard_parallel; standard_parallel1; standard_parallel2;
longitude_of_central_meridian; longitude_of_projection_origin; latitude_of_projection_origin; > straight_vertical_longitude_from_pole; longitude_of_prime_meridian; semi_major_axis;
inverse_flattening; earth_radius; +proj; +x_0; +y_0; +k_0; +k_0; +lat_1; +lat_1; +lat_2;
+lon_0; +lon_0; +lat_0; +lon_0; +lon_0; +a; +rf; +a

> rrqpe
class       : RasterLayer
dimensions  : 5424, 5424, 29419776  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -0.5, 5423.5, -0.5, 5423.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc
names       : ABI.L2..Rainfall.Rate...Quantitative.Precipitation.Estimate
zvar        : RRQPE

Then I set the CRS:
> rrqpe@crs <- sp::CRS('+proj=geos +h=35786023.0 +a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs')
> rrqpe
class       : RasterLayer
dimensions  : 5424, 5424, 29419776  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -0.5, 5423.5, -0.5, 5423.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=geos +h=35786023.0 +a=6378137.0 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs
data source : OR_ABI-L2-RRQPEF-M6_G16_s20202340210217_e20202340219525_c20202340220028.nc
names       : ABI.L2..Rainfall.Rate...Quantitative.Precipitation.Estimate
zvar        : RRQPE

Then I try to warp (reproject) the image to EPSG:4326
> rrqpe_wgs84 <- projectRaster(rrqpe, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84")

But something is wrong. The extent is tiny. I was expecting something like half world of data.
> rrqpe_wgs84
class       : RasterLayer
dimensions  : 5436, 5436, 29550096  (nrow, ncol, ncell)
resolution  : 8.98e-06, 9.04e-06  (x, y)
extent      : -89.50005, -89.45123, -4.774882e-05, 0.04909369  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /tmp/RtmpfHzSHu/raster/r_tmp_2020-08-21_231812_20555_52886.grd
names       : ABI.L2..Rainfall.Rate...Quantitative.Precipitation.Estimate
values      : -49.68673, 49.87008  (min, max)

I want to plot the data (rrqpe_wgs84) on a world map like world["continent"]["geom"] from library spData. Can someone help me on this reprojection using R libraries? Although I'm not using here, I have rgdal installed. I'm trying to avoid using many file operations with gdalUtils on R 3.4.4.


texts: Are there any GIS component that can be used for developing GIS applications for Android phones?

texts: I know that ESRI has ArcGIS APIs for Android what other options are there for me as developer, if I'm about to develop GIS applications for an Android phone?


texts: GRASS module g.mremove returns irrational outputs

texts:
I want to use g.mremove to remove multiple files in my mapset,

import os
import sys
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

g.run_command('g.mremove', rast = 'b_*')
print g.read_command('g.list', _type='rast')

but the output still has b-beginning files,

*** Remote Interpreter Reinitialized  ***
>>>
----------------------------------------------
raster files available in mapset <TC>:
acc        b16        b23        b30        b7         bas3       basin_test
b10        b17        b24        b31        b8         bas4       dra
b11        b18        b25        b32        b9         bas5       hbas
b12        b19        b26        b33        bas        bas6       hc
b13        b20        b27        b4         bas1       bas7       str
b14        b21        b28        b5         bas10      bas8       vis
b15        b22        b29        b6         bas2       bas9

raster files available in mapset <PERMANENT>:
hc

----------------------------------------------

How to fix this problem？ I am using GRASS 6.4.3, python 2.7.6 under Windows 8.1 x64.

Update#1
I have edited the g.run_command line above to this,

g.run_command('g.mremove', flags = 'f', rast = 'b_*')

the output is the same, then I try to do it in the python shell of the GRASS layer manager, and here are the code and output,

>>>g.run_command('g.mremove', flags = 'f', rast = 'b_*')
0
>>>g.read_command('g.list', _type = 'rast')
----------------------------------------------
raster files available in mapset <TC>:
acc        b16        b23        b30        b8         bas5       hc
b1         b17        b24        b31        b9         bas6       str
b10        b18        b25        b32        bas        bas7       vis
b11        b19        b26        b33        bas1       bas8
b12        b2         b27        b4         bas10      bas9
b13        b20        b28        b5         bas2       basin_test
b14        b21        b29        b6         bas3       dra
b15        b22        b3         b7         bas4       hbas

raster files available in mapset <PERMANENT>:
hc

Still give me irrational outputs.


texts: Error in qgis:selectbylocation python processing method

texts:
I am having trouble with my code for selecting all points within a buffer layer, with the error given below. I am using this documentation as reference for the field formats:
https://docs.qgis.org/2.6/en/docs/user_manual/processing_algs/qgis/vector_selection_tools/selectbylocation.html . I also tried the solution posted here How to use qgis:selectbylocation in PyQGIS? . Is there something that I'm missing?

processing.runalg("qgis:selectbylocation",
                  "/Users/pgcseismolab/Desktop/qgis/moment_all.shp",
                  "/Users/pgcseismolab/Desktop/jaden_pgc_archive/Blast_Files/Blast_buffers/blast_induced_buff.shp",
                  False,
                  False,
                  False,
                  0)

Error: Wrong number of parameters
ALGORITHM: Select by location
    INPUT <ParameterVector>
    INTERSECT <ParameterVector>
    PREDICATE <ParameterGeometryPredicate>
    PRECISION <ParameterNumber>
    METHOD <ParameterSelection>
    OUTPUT <OutputVector>

METHOD(Modify current selection by)
    0 - creating new selection
    1 - adding to current selection
    2 - removing from current selection

EDIT: I am using QGIS 2.14, and have updated my code as follows (with the error below). Following documentation posted here: http://docs.qgis.org/2.14/en/docs/user_manual/processing_algs/qgis/vector_selection_tools.html?highlight=selectbylocation

processing.runalg("qgis:selectbylocation",
                  "/Users/pgcseismolab/Desktop/qgis/moment_all.shp",
                  "/Users/pgcseismolab/Desktop/jaden_pgc_archive/Blast_Files/Blast_buffers/blast_induced_buff.shp",
                  'within',
                  0)

Error: Wrong number of parameters
ALGORITHM: Select by location
    INPUT <ParameterVector>
    INTERSECT <ParameterVector>
    PREDICATE <ParameterGeometryPredicate>
    PRECISION <ParameterNumber>
    METHOD <ParameterSelection>
    OUTPUT <OutputVector>


texts:
In one of my work, I just need the QA band of Landsat, but USGS Earth Explorer does not provide band by band data download opportunity. As I am interested to get all the available Landsat QA data from Landsat 1 through 8, downloading full scene(all bands) will be useless for me.

I have tried SCP plugin with QGIS, did not found any solution.

Is there any way I could download only QA band and the metadata rather than downloading the whole scene?


texts: Reasons to use plain synchronous .getInfo() in GEE

texts:
When querying information from EE servers it is advised to use asynchronous callbacks, e.g. .evaluate() or .getInfo(callback). So the question is: are there any reasons to use plain synchronous .getInfo() at all, aside from its siplicity?

P.S. Of course there are cases where some parts of code must be executed after retrieving information from server, but even in those cases I see a better option in putting that code in a callback function.


texts: How do I calculate polygon area with arcobjects?

texts:
I'd like to develop an application with ArcEngine that can calculate the area of a polygon.

I wrote the below code to select a polygon in a combobox:

for (int i = 0; i < pmapcontrol2.LayerCount; i++)
{

    Pfeaturelayer = (IFeatureLayer)pmapcontrol2.get_Layer(i);
    if (Pfeaturelayer.FeatureClass.ShapeType == esriGeometryType.esriGeometryPolygon)
    {

        cmbfeaturelayer.Items.Add(pmapcontrol2.get_Layer(i).Name);
    }
}

But I couldn't write code to calculate pfeaturelayer area.
I think I should use Ilayer interface and Ilayer use Ifeature to calculate area not ifeaturelayer.

Is there any way to use Ifeature for this application?


texts:
In the Browser it's not possible to delete Shapefiles. In ArcMap you can delete Shapefiles easily within the Catalog.

How can I delete Shapefiles in QGIS without minimizing the QGIS-Window, navigating to the folder on my computer, select all the files that belong to a Shapefile (shp, dbf, shx, prj, etc.) and delete it?


texts:
When I add a new layer to a QGIS project it is often not visible in the viewing pane, although it is visible in the side Layers Panel.

I would like to make the layer visible and have tried to use the mouse to drag and drop the layer to the top of the layers panel.  This however crashes the programme.  The only work around I have found is to remove the basemap raster layer and then reinstall it.

I have attached a screenshot which shows the 'boundary' layer in the layers panel but not visible in the main viewing pane.


texts: What does it mean to reproject a satellite image from utm (zone 13N) to wgs84?

texts:
I'm trying to understand some legacy code that processes landsat satellite images. Source images are in UTM and they are reprojected to WGS84. I'm a GIS noob but I understand that WGS84 is a datum for a geographical coordinate system and UTM is a projected coordinate system so what does it mean to reproject from UTM to WGS84? Does it mean to reproject to some standard projection of WGS84? What is that standard?

Please see following code:

# First open raw landsat SR image
source_file = os.path.join(VALHALLA_ROOT, 'test/data/2015_035024/LC08_L1TP_035024_20150626_20170226_01_T1_sr_band1.tif')

with rasterio.open(source_file) as src:
    data1 = src.read()
    print(data1.shape)
    print(src.meta)

  (1, 7901, 7801)

  {'count': 1, 'driver': 'GTiff', 'height': 7901, 'width': 7801, 'transform': Affine(30.0, 0.0, 518985.0,
         0.0, -30.0, 5846715.0), 'dtype': 'int16', 'nodata': -9999.0, 'crs': CRS({'init': 'epsg:32613'})}

# Now reproject to WGS 84
from rasterio.warp import calculate_default_transform, reproject, Resampling

def reproject_tif(source_file, destination_file):
    """Re-projects tif at source file to WGS84 at destination file.

    Args:
        source_file: file to re-project
        destination_file: file to store re-projection

    Returns:
        destination_file: where the re-projected file is saved at
    """

    if os.path.isfile(destination_file):
        logging.info("Projected file already exists, skipping reprojection")
        return

    with rasterio.open(source_file) as src:
        dst_crs = 'EPSG:4326'
        transform, width, height = calculate_default_transform(
            src.crs,
            dst_crs,
            src.width,
            src.height,
            *src.bounds
        )

        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })

        with rasterio.open(destination_file, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest,
                    num_threads=-1
                )

        return destination_file

dest_file = os.path.join(VALHALLA_ROOT, 'test/data/2015_035024/Reproject_LC08_L1TP_035024_20150626_20170226_01_T1_sr_band1.tif')

with rasterio.open(dest_file) as src:
    data2 = src.read()
    print(data2.shape)
    print(src.meta)

  (1, 6140, 9769)

  {'count': 1, 'driver': 'GTiff', 'height': 6140, 'width': 9769, 'transform': Affine(0.00035596559131628866, 0.0, -104.73152062773751,
          0.0, -0.00035596559131628866, 52.7699321953733), 'dtype': 'int16', 'nodata': -9999.0, 'crs': CRS({'init': 'epsg:4326'})}

%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(2, figsize=(10, 9))
ax[0].imshow(data1[0, :, :], cmap='gray')
ax[1].imshow(data2[0, :, :], cmap='gray')


texts: Google Earth graphics card problem in Ubuntu 14.04 64bit

texts:
I am struggling to get Google Earth working properly in Ubuntu 14.04 64bit. The program opens correctly and does not crash, as many report over the web. But it has trouble working with the graphics card, returning a "clamped polygons not supported" message.

(The map in question was rendered properly in Windows in a dual-boot machine, so the card itself is not the problem, it would seem.)

I downloaded/installed the 64bit .deb and also the 32-bit libraries that are required, as instructed here and also here. I also tried installing the 32-bit GE with even less success.

Has anyone solved this issue? I am now inclined to install the win version of GE over wine, but wish to ask for advise first.

Or might installing an earlier GE version do the trick, as suggested here?

Some possibly relevant info about my graphics card follows

$ lspci -vnn | grep VGA -A 12
00:02.0 VGA compatible controller [0300]: Intel Corporation 2nd Generation Core Processor Family Integrated Graphics Controller [8086:0116] (rev 09) (prog-if 00 [VGA controller])
Subsystem: Dell Device [1028:052e]
Flags: bus master, fast devsel, latency 0, IRQ 46
Memory at f0000000 (64-bit, non-prefetchable) [size=4M]
Memory at e0000000 (64-bit, prefetchable) [size=256M]
I/O ports at 2000 [size=64]
Expansion ROM at <unassigned> [disabled]
Capabilities: <access denied>
Kernel driver in use: i915

00:16.0 Communication controller [0780]: Intel Corporation 6 Series/C200 Series Chipset Family MEI Controller #1 [8086:1c3a] (rev 04)
Subsystem: Dell Device [1028:052e]
Flags: bus master, fast devsel, latency 0, IRQ 47


texts:
Topo to Raster conversion tool in the 3D Analyst toolset, It gives me the following error and tells me to change the file name, which doesn't work
ERROR 010168: Output C:\WINDOWS\system32\zpa09 already exists. Failed to execute (TopoToRaster).


texts: System variable PATH overwritten in QGIS?

texts:
The system env variable 'PATH' seems to be overwritten in Qgis, replacing it with their own values. Not sure if it was intended or not.

Anyway in version 2.2, I could solve it by go into settings and appending the PATH variable with the values copied from system properties => environment variables.

Apparently for version 2.4, the above solution does not work anymore and Qgis will always overwrite the PATH variable with their own values.

I needed the PATH variable as my project requires the execution of external programs from Qgis by making use of python subprocess module. Now in v2.4, external programs cannot be executed. I do not wish to define absolute path for all the external programs as the project will frequently be ported to different systems.


texts: ArcObjects - Create and Draw a Transparent Polygon

texts:
I'm trying to create a transparent Polygon in ArcMap; however, it's never transparent.  Here's the code I have.

// Create polygon
IPolygon polygon = new PolygonClass();
// create & add points,......

// Create lineSymbol
ISimpleLineSymbol lineSymbol = new SimpleLineSymbolClass();
// set color & style,....

// Create fillSymbol
ISimpleFillSymbol fillSymbol = new SimpleFillSymbolClass();
fillSymbol.Outline = lineSymbol;
fillSymbol.Color = //set to color red//;
fillSymbol.Color.Transparent = 0;
fillSymbol.Style = esriSimpleFillStyle.esriSFSSolid;

// Draw Shape
myMap.DrawShape(polygon, fillSymbol);

RESULT ==> A solid red polygon is drawn.  The Transparent setting is ignored.

Thoughts?


texts: How to not process when query is outside the raster extent ? (Error when Completely Outside)

texts:
Given a West North East South bounding box, I want to crop this data from a hundreds of inputXX.tif. Nearly every times but once, the query first print out that it "falls completely outside raster extent but is going on [processing it] however", then uselessly start to process the out-of-range raster, spenting huge amoint of time on EACH file.

$gdal_translate -projwin 0 55 -10 45 ./data/input${i}.tif ./tmp/crop.${i}.tmp.tif
Input file size is 6001, 6001
Computed -srcwin 296400 -57000 38400 39000 from projected window.
Warning: Computed -srcwin 296400 -57000 38400 39000 falls completely outside raster extent. Going on however.
0...10...20...30...40...50...60...70...80...90...100 - done.

As I'am processing SRTM3's 871 files via a makefile-bash loop...

#!bash
files=(`ls ../data/cgiar-csi_90/ | grep .tif$`)
for i in "${!files[@]}"
do
    gdal_translate -projwin $WEST $NORTH $EAST $SOUTH ../data/cgiar-csi_90/${files[${i}]} ./tmp/crop.${i}.tmp.tif
done
gdal_merge.py -o ./crop_xl.tmp.tif ./tmp/crop.*.tmp.tif -n -32768

...How to NOT process such raster when query is outside the raster extent ? and pass without error.

Note: gdal_translate doesn t seems to have native way to jump over such out-of-range inputs, as far as i understand.


texts:
I have read documentation on overwriting ArcGIS service via python. The overall process to programmatically overwrite service is to save a service definition draft, in the .sddraft file the set service type to from esriServiceDefinitionType_New to esriServiceDefinitionType_Replacement, then stage the .sddraft to a .sd file, and upload the .sd file.

What is the difference between setting the service type to from esriServiceDefinitionType_New to esriServiceDefinitionType_Replacement in the .sddraft file and setting it in the .sd file?

I ask this because I have several .sd files already created, and with a simple 7zip extraction I can open the .sd file and see the service type setting in the manifest.xml.

In theory I would like to make the service type change here in the .sd file, then upload and overwrite my existing service.

My current plan is to continue testing, and if all seems well, stick to this approach.

My environments: ArcGIS Server 10.3.1 & ArcMap 10.3.1


texts:
I have 6 samples of water quality parameters taken on a 10 km stretch downstream (Cooum River in Chennai), so it's basically the estuary.

The river discharge and flow is extremely low, therefore I assumed the current to be not an important factor.

Which interpolation method would make sense in this set up?

Since the number of sample points is extremely low I moved away from IDW.

At the moment I am considering Indicator Kriging, which was recommended for pollution data in general.


texts: How to get pixel values/cell values for stored raster files in postgis?

texts:
Raster files in tiff format are storing into database.However I gets only an attribute table.
Is there any way to retrieve pixel values/cell values in the table.

I have used raster2pgsql command for storing.


texts: Fixing Topology (Gaps and Overlaps) in QGIS

texts:
I am digitizing a new polygon shapefile of property boundaries based on orthoimagery using QGIS Zanzibar. I'm acquainted with topology editing in ArcMap using geodatabases. I did so with snapping to vertices and segments enabled.

Using the Topology Checker Plugin in QGIS 3.8 Zanzibar and in 3.6 Noosa I can identify digitizing errors such as gaps and overlaps.

However when attempting to use v.clean in the Processing window I keep running into the following errors:

WARNING: Concurrent mapset locking is not supported on Windows
WARNING: Column 'id' : type int8 (bigint) is stored as integer (4 bytes) some data may be damaged
WARNING: Column 'id' : type int8 (bigint) is stored as integer (4 bytes) some data may be damaged
WARNING: Unable to create index
WARNING: Unable to copy table <Block_6_Parcels_Test> for layer 1 from <vector_5d9f8b040caf72@PERMANENT> to <output39ece58e87d34cfa8e0b689616f7aafc>
WARNING: Failed to copy attribute table to output vector map
WARNING: Number of centroids exceeds number of areas: 7 > 5
WARNING: Number of incorrect boundaries: 2
WARNING: Number of centroids outside area: 2
ERROR: Key column <> not found
WARNING: No areas found, but requested to be exported. Will skip this feature type.
WARNING: No attribute table found -> using only category numbers as attributes
WARNING: Output layer is empty, no features written

I have two fields that could be used as key fields, either the id column that is the default generated first column when you create a new shapefile, or the PARCEL_ID. Both are created as @row_number in field calculator. Both are populated.

I've tried just doing one error rule at a time (i.e. fixing gaps first with v.clean only, then fixing overlaps with v.clean) but both fail. I've tested this in both 3.6 Noosa and 3.8 Zanzibar.

How should I go about fixing my topology? A link to the shapefile https://www117.zippyshare.com/v/vhZODckl/file.html


texts: I have all the settings right to generate pictures using Atlas, but I can't access the 'Controlled by Atlas' checkmark in the print layout. It's greyed out, and whatever I try I'm unable to change this. I have a shapefile as layer, and an XLS file as Coverage layer. No background. Any idea what I'm doing wrong?


texts: User memory limit exceeded error for classification in Google Earth Engine

texts:
I am trying to do classification using cart. Given below is the code that was used, but I am getting 'user memory limit exceeded' error. Is there no other way to rectify this other than by reducing the extent of the ROI, randomly sampling or by knocking out a few bands? Am I doing something wrong here? This is the asset that I am using to train: https://code.earthengine.google.com/?asset=users/RohitNandakumar/train

var sent = ee.ImageCollection("COPERNICUS/S2_SR"),
    table = ee.FeatureCollection("users/RohitNandakumar/train"),
    geometry =
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[79.33418599339291, 29.71454696770813],
          [79.33418599339291, 29.317189386242045],
          [79.86496297093197, 29.317189386242045],
          [79.86496297093197, 29.71454696770813]]], null, false);

var land = sent.filterDate('2019-04-1','2019-11-01')
.filterMetadata("MGRS_TILE","equals",'44RLT')
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE",5))
.sort('system:time_start');
var image = ee.Image(land.toList(land.size()).get(1)).clip(geometry);
var label = 'MC_ID';
var bands = ['B2','B3','B4','B5','B6','B7','B8','B8A','B9','B11','B12']
var training = image.select(bands).sampleRegions({
  collection: table,
  properties: [label],
  scale: 10
});
var trained = ee.Classifier.cart().train(training, label, bands);
var classified = image.select(bands).classify(trained);
Map.addLayer(classified,false);


texts:
Depending on this example

I've created the following WPS request to get the nearst airport from a given point.

<?xml version="1.0" encoding="UTF-8"?>
<wps:Execute version="1.0.0" service="WPS" xmlns:tst="tst" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>gs:Nearest</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>features</ows:Identifier>
      <wps:Reference mimeType="text/xml; subtype=wfs-collection/1.0" xlink:href="http://geoserver/wfs" method="POST">
        <wps:Body>
          <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
            <wfs:Query typeName="tst:ne_10m_airports"/>
                 <ogc:Filter>
                     <ogc:PropertyIsGreaterThan>
                          <ogc:PropertyName>scalerank</ogc:PropertyName>
                          <ogc:Literal>5</ogc:Literal>
                     </ogc:PropertyIsGreaterThan>
                 </ogc:Filter>
          </wfs:GetFeature>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>point</ows:Identifier>
      <wps:Data>
        <wps:ComplexData mimeType="application/json">
        <![CDATA[
        {
          "type": "Point",
          "coordinates": [
              21.5,
              15.5
          ]
        }
        ]]>
      </wps:ComplexData>
      </wps:Data>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>crs</ows:Identifier>
      <wps:Data>
        <wps:LiteralData>EPSG:4326</wps:LiteralData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="application/json">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

The filter inside the WFS:Query is not executed, although I put a criteria scalerank > 5, the WPS request returns an airport with scalerank = 4 !!

What am I missing or doing wrong?


texts: Resolution and other characteristics of multipatch

texts:
i'm trying to understand how multipatch works. I have one element and i would like to change its resolution and position. when i try to change the position, there are limits for that, only can move this in Y 10 units ande X-Z 0.1 units ¿someone know how can i move it, are the limits adjustable?. i'm thinking it's probably that the domain has a important roll in this case, but i couldn't modify this one ¿is that possible?

When i try to change the resolution, in the moment of import 3D files i can to modify XYZ resolution and tolerance, but then when i check this items in the multipatch, they dont have the values introduced by me. ¿are these characteristics adjustable?

i don't know, perhaps these characteristics are stipulated by default depending on kind of archive and they are not adjustable. ¿Does anyone have any ideas about this?


texts: Increasing raster exent with Raster Calculator is modifying band range

texts:
I have a TIFF with a single band 16bit signed. I used the Raster Calculator tool to extend its extent:

Con(IsNull("mytif.tif"), 1, "mytif.tif")

Using ArcGIS Pro 2.0.0 I set the environment to the extent "As Specified Below" and the Snap Raster to the original TIFF (so that the cells align) and set the cell size to that of the original TIFF as well.

The original TIFF band has a range from 1 to 2487, the output raster from the tool has -132 to 2684.

Why are the values changing and how can I change them back to their correct scale (the represent elevation data at a specific resolution)?


texts: How to select centroid of multiple nearby points with same attribute

texts:
This question is very similar to #131823, that is not solved yet.

I want to merge several points like bus stops into one single point in the centroid of all selected points. The bus stops have a attribute "name". The problem is that these names are not necessarily unique worldwide, so there needs to be a spatial limit if several points can be combined to a MultiPoint. This done, we can calculate the centroid with PostGIS: http://www.postgis.org/docs/ST_Centroid.html

How can we solve this problem with Postgis functions?


texts: Comparing and reporting shift between two shapefiles using ArcGIS Desktop?

texts:
I have two shapefiles.  One is an original route polyline file and another is an updated polyline file, where the shape is modified to google maps.

How can I compare the updated shapefile with the original and report how much shift or variation is there in the second file, with reference to the original ?

There seems to be a feature compare tool in ArcGIS , but seems that it will not suit the purpose.


texts: Removing tooltip label border completely in Leaflet.js map?

texts:
I am trying to completely remove a label boundary/border (bubble) that by default appears when using my labels of a point GeoJSON layer.

The layer I have is defined as below:

var points_layer = new L.geoJson(points, {
    onEachFeature: onEachPointFeature,
    pointToLayer: function (feature, latlng) {
         return L.circleMarker(latlng, labelMarkerOptions);
    }
});

the options for the label are

var labelMarkerOptions = {
  opacity: 0,
  fillOpacity: 0,
  fillColor:'black',
  fill: 'false',
  interactive: false

};

the oneachFeature function  is

function onEachPointFeature(feature, layer) {
  if (feature.properties) {
     layer.bindTooltip(feature.properties.v2_number.toString(), {permanent: true, className: 'labelstyle'});
   }
}

using the style below I was able to remove almost all of the borders around the label text but a small white triangle remains next to the label which I cannot get rid of. This somehow seems to be a svg path or something (but can determine what it really is fully or where it comes from exactly).

How can I accomplish this?

.labelstyle {
      all: revert;
    color: green;
    font-size:14px;
    font-weight: 700;
    fillColor: none;
    fillOpacity: 0;
    background-color: none;
    border-color: none;
    background: none;
    border: none;
    box-shadow: none;
    margin: 0px;
    cursor: none;
    direction: 'center';
    interactive: false;
    fill: false;
}

Below is how the remaining border of the original bubble is currently looking. The White triangles to the right of the green labels is (what is there by default) and what I would like to get rid of

I do not think that this is a duplicate of Overriding Leaflet tooltip style?
 because that is where I started using the tips in there but still could not get the white triangle styled at all.

Anyone has insight how to alter my css to make the triangle disappear ?


texts: Infinity Downslope Influence with TauDEM

texts:
Does anybody now how to calculate Infinity Downslope Influence with the TauDEM version?

the command was present until version 3.x but now disappeared I don't know if it was renamed or simply removed.


texts: Using a ArcGIS Server Feature Service's addFeature method: Parameters?

texts:
I am wanting to post data back to a feature service using the addFeatures method on a feature service.  Here it says to use a query and an esri rep says to use a query task.  But my "Query" object doesn't seem to have a "features" definition as suggested in the linked doc page above where they say the parameters table.  Anyone have any ideas on how to do the addFeatures.

QueryTask queryTask =
                new QueryTask("http://server/ArcGIS/rest/services/geoprocessing/featureServiceName/FeatureServer/0/addFeatures");

            ESRI.ArcGIS.Client.Tasks.Query query = new ESRI.ArcGIS.Client.Tasks.Query();
            query.???


texts: What graph traversal algorithm should I use?

texts:
I would like to write an algorithm which can traverse a graph and hopefully later I can implement it to use for an indoor navigation system. The graph would come from floor plans of a building and the graph nodes represent the building objects such as doors, corridors, stairs, rooms, etc.

The navigation system will accept different users with different capabilities, these capabilities are stored in user profile. For each individual user (or a group of users of the same type, e.g. staff member) the system returns a path specific to that user/s considering his personal functionality. For example if the user of the system is using a wheelchair the system should exclude all the stair nodes from the path (they will be marked non-traversable). Similarly if in order to enter into an indoor space, the user needs to have a key or be an staff member but he is NOT, the system should mark the node (e.g. the door which requires the key) non-traversable and should search for those paths which don't require these specifics requirements.

This was the first part and here I use a function to check these conditions for a given user. The inputs to this function are a node and a user. If the conditions that are required by the node are NOT fulfilled (means user profile doesn't match the condition required by the node) the function returns false= non-traversable. If the condition holds true, the node stays isTraversable and will be stored somewhere and will be considered for the final path calculation.

The second part is that I have a list that contains all the nodes that are traversable and now the shortest path will be calculate based on 3 different costs, Energy, time and money. The final cost can be a combination of these costs and it sets the weights for the edges.

My problem is with the first part. I want to know how can I make an algorithm that produces the subgraph (the list with the traversable nodes in it). And later how can I use this algorithm to modify Dijkstra to fit this system.

Can anybody help with writing a pseudocode for the algorithm so that I can understand how would it work in clear steps or any hints on how can I make these steps into pseudocode?

Please note that at this stage the performance is not very important for me.

Thanks in advance!


texts: Is there a way to apply an arbitrary coordinate translation to a layer in QGis 2.8.2? I have some data sets that are in a truncated coordinate system, and I need to shift them by a certain amount to line up with the real world coordinate system.


texts: ArcGIS Issue with Python wildcard in Label Expression?

texts:
Here is my current expression:

def FindLabel ( [FULL_NAME] ):
fn = [FULL_NAME]
if "COUNTY ROAD" in fn:
 return  "CR " + fn[11:]
elif "% COUNTY ROAD" in fn:
 return fn[:1] + "CR" + fn[13:]
else:
 return fn

This is a valid expression but the line - elif "% COUNTY ROAD" in fn:
     return fn[:1] + "CR" + fn[13:] isn't working. Some county roads have a pre directional (i.e. E or W) then COUNTY ROAD XXXX. I am trying to return that pre-directional with % wildcard if it exists in the FULL_NAME field.

Is this not correct?


texts: Possible to show aggregate votes in a tooltip on a markercluster?

texts:
I'm building an interface that would take poll data in and show some results over a marker cluster using Mapbox but I don't have any experience with it.

The interface would look much like https://www.mapbox.com/mapbox.js/example/v1.0.0/leaflet-markercluster/

If you click on the cluster instead of zooming in a tooltip would show up with a breakdown of (yes/no) votes in a pie chart like this https://www.dropbox.com/s/hv4c8708wlhe4mt/Screenshot%202014-11-30%2015.27.53.png?dl=0

In the GeoJSON data I was thinking each data point would have a property to state the vote value and somehow I could look at them in aggregate at the cluster level

geoJSON structure

{
  "type": "Feature",
  "geometry": {
    "type": "Point",
    "coordinates": [125.6, 10.1]
  },
  "properties": {
    "vote": "yes" //no
  }
}

Is this possible using these libraries and geoJSON data structure? Am I on the right track?


texts: COGO rotation skewed, likely coordinate issue?

texts:
I am COGO'ing a survey plan from June of 1972, the scan of the paper plan is not so great and does not list any Datum or declination used, but only a scale bar showing measurements are in feet, and the coordinates are DD.MM.SS

I entered a ground correction to get my feet into meters, and when COGO'ing it I wind up with a shape which closely aligns to what the cadaster says it should look like, but is obviously rotated SW.  I suspected this is to do with some coordinate system mismatch so wondering for some advice.

I believe in 1972 there was a new variant of WGS that came out, so not sure if the plan might be in either WGS1966 or WGS1972, but can't even guarantee it is in WGS, only making a guess based on the DD.MM.SS coordinates, but for all I know this is using NAD27...

Below is a picture of the COGO (green line) and where it sits with the cadaster, the POC is the NW corner and I had to COGO two lines in each direction as the area where they would close the traverse is indecipherable on the paper scan.

To complicate things, the POC for the traverse on the paper scan is very hard to read, but putting in the rough coordinates I can read doesn't put me near the shape in question, instead it puts me several kilometers away in the middle of water, not on land, but I figure the POC could be a key to sorting out the Datum for the rest of it and if I can figure out how to make those coordinates fall somewhere close to the COGO'd shape I'll be on the right track.  Here's a pic of the POC:

Any advice on how to correct for this unknown coordinate system and make this COGO shape line up?

Should I be adjusting all bearings with declination from june 1972?


texts:
I am looking for a simple point shapefile which indicates where the cities of the world are.
I cannot seem to find one which is free of charge and meets my requirements.
Ive got this one: http://techcenter.jefferson.kctcs.edu/data/ but I want smaller cities too...

Any ideas where I can get them?


texts: Programmatically stop and start ArcGIS Server map services using Python in ArcGIS 10.5?

texts:
In his answer to Starting and stopping ArcGIS Server map services using Python? @Jeff Berry posted how can we can programmatically stop and start ArcGIS Server 10.0 map services using Python.

How do I stop services in ArcGIS 10.5?

I couldn't find AGSSOM there.


texts:
I would like to provide the user of a plugin with progress status using messages and a progress bar.

My plugin processes several files using memory files, processing, writing to disk etc.  It all works well, but takes 10 to 90 seconds run.  In QGIS 2.18 it ran for about 5 seconds displaying progress text in the status bar, then the Windows wait icon displays and QGIS is unresponsive until complete - status bar hangs and QGIS looks frozen.  In version 3.0, (the same plugin after updating the code), QGIS appears to freezes immediately (although the code is working properly), and status text does not display at all until after the routine has finished.

I am using Windows 10 and QGIS 3.4.2.  My code takes long enough to require a status bar, but my attempts to display any information only occurs after the code has finished.  I have looked into tasks and running code in the background, but there appear to be bugs yet to be resolved  (QgsTask and task manager cause crash when used in function scope).  Is there any way in QGIS 3.0 to show progress during the processing steps?

Update: I have had some success with a progress bar as a widget in the message bar.  By adding self.iface.mainWindow().repaint() after setting the bar value (setValue) the bar will display properly.  If the code takes more than 5 seconds the screen freezes and no updates to QGIS until the code has run.  If I add time.sleep(0.1) before and after setting value and repainting, the status bar works properly.

Is it normal to have to repaint the screen every time the status bar is to increment?  Is the only workaround to "jiggle" the screen with the sleep-timer so it won't freeze?

The above solution came from Unresponsive MessageBar in QGIS 3


texts: Sorting out curved and straight lines in QGIS

texts: I have a layer in QGIS with thousands of lines and I would need to sort out the curved lines and lines that are straight. Can I do this easily in QGIS?


texts: arcpy SelectLayerByAttribute_management runs in the Python window in ArcGIS Pro, but fails to run from Spyder or PyCharm

texts:
I have this script for querying a feature layer and then copying to a geodatabase. It works when I run it in the ArcGIS Pro Python window, but not when I run it from Spyder or PyCharm. When I run in Spyder the script fails at SelectLayerByAttribute_Management. I have tried multiple variations on the structure of the definition query. Hopefully someone can advise :)

The field the where_clause is applied to is a text field.

The error I am getting from Spyder is:

  runfile('C:/Users/LENOVO-580/REMI/Join_Data.py',
  wdir='C:/Users/LENOVO-580/REMI') Made feature layer Traceback (most
  recent call last):

  File "", line 1, in
      runfile('C:/Users/LENOVO-580/REMI/Join_Data.py', wdir='C:/Users/LENOVO-580/REMI')

  File "C:\Program
  Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\spyder_kernels\customize\spydercustomize.py",
  line 704, in runfile
      execfile(filename, namespace)

  File "C:\Program
  Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\spyder_kernels\customize\spydercustomize.py",
  line 108, in execfile
      exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/LENOVO-580/REMI/Join_Data.py", line 61, in
      Selection = arcpy.SelectLayerByAttribute_management(out_layer, "NEW_SELECTION", where_clause)

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\management.py", line 7238, in
  SelectLayerByAttribute
      raise e

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\management.py", line 7235, in
  SelectLayerByAttribute
      retval = convertArcObjectToPythonObject(gp.SelectLayerByAttribute_management(*gp_fixargs((in_layer_or_view,
  selection_type, where_clause, invert_where_clause), True)))

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing_base.py", line
  496, in
      return lambda *args: val(*gp_fixargs(args, True))

  ExecuteError: ERROR 000358: Invalid expression Failed to execute
  (SelectLayerByAttribute).

import arcpy

monthYear = 'October 2018'

#feature from path
in_features = r'Z:\Data_Path{0}\Map
Data\GDB_{0}.gdb\Assessed_Areas'.format(monthYear)

# Set the workspace
arcpy.env.workspace = r'Z:\Data_Path{0}\Map Data\GDB_{0}.gdb
out_layer = "out_layer"
where_clause = "Relevant_Column = '0'"
output = r'Z:\Data_Path{0}\Map Data\GDB_{0}.gdb\Output'.format(monthYear)

arcpy.MakeFeatureLayer_management(in_features, out_layer)
print("Made feature layer")
Selection = arcpy.SelectLayerByAttribute_management(out_layer,
"NEW_SELECTION", where_clause)
arcpy.CopyFeatures_management(Selection, output, None, None, None, None)


texts: Osmfilter - filter every node with at least one tag

texts:
I'm using osm data and the osmfilter tool but I can't achieve to filter all nodes with at least one tag. I want to drop nodes with no tags.

My approch is:

osmfilter DataOSM.o5m --keep-nodes= --keep-tags="*=" -o=test.osm

But it does not work. All nodes are kept.


texts:
Within the openlayers documentation I’ve noticed that “options” is many times used to show that the constructor takes parameters. But where can I get the details for this parameter {object}? How does one know what parameters the constructors accepts?

For example: below I have put the constructor documentation for OpenLayers.Filter.Spatial.

How do I know what properties I can put within {Object}?


texts: GDAL C++ Re-Project Coordinate with Elevation (Z)

texts:
I need some assistance in determining why the following code is not re-projecting the elevation parameter. It works perfectly for the Latitude and Longitude portions but seems to completely ignore the Z parameter that's being passed in.

GDALAllRegister();

OGRSpatialReference srFrom;
srFrom.SetProjCS("WGS 84 / UTM zone 33N");
srFrom.SetWellKnownGeogCS("EPSG:32633");
srFrom.SetUTM(33,TRUE);
auto srFrome = CPLGetLastErrorMsg();

if (!std::string(srFrome).empty())
{
    auto p = "";
}

OGRSpatialReference srTo;
srTo.SetWellKnownGeogCS("EPSG:4326");
auto srToe = CPLGetLastErrorMsg();

OGRCoordinateTransformation* coordTrans =
OGRCreateCoordinateTransformation(&srFrom, &srTo);
auto cte = CPLGetLastErrorMsg();

double x1 = 1291784.057793292;
double y1 = 2724159.114436987;
double z1 = 453.896;

int reprojected = coordTrans->Transform(1,&x1,&y1,&z1);

//Z always returns the original value passed in.


texts: Is it possible to open a shapefile using a ArcGIS Runtime for Windows Mobile application

texts: I need to build a tablet application using ArcGIS Runtime for Windows Mobile. Is there a way to open a shapefile (or a pgdb) without having to convert it to a .db file first?


texts: QGIS Select By Location geometric predicates explanations?

texts:
Im trying to test if a polygon either a) touches or b) is only semi-contained by another polygon. The first is self-explanatory i guess, but the second could be many of the different predicates offered in QGIS.

I was therefore wondering if there are any more indepth explanations of the geometric predicates used in QGIS' select by location. (the only slightly basic explanation is Documentation "Join Attributes by Location" QGIS 2.18 which doesn't give any real clue as to what the terms cover)


texts: How to configure GeoExplorer to render as HTML the content of the "Table panel" cells (especially as URL)?

texts:
I want to show in GeoExplorer the values of feature attribute set as url link, e.g.
 link
to give users the chance for clicking them.
It works inside the feature info window, but in the table panel (which is opened after clicking the query button) the same attribute value is displayed as plain text.
How can i change it?

Regards
Simon


texts: Is extrusion of a line (which is curved) into a plane at a certain angle possible in ArcGIS?

texts:
I have a polyline shape file which is a fault trace on earth's surface.

I want to extrude this into a 3-D plane at a certain angle (which is dip angle of the fault).

However, the extruded plane should start somewhere below the ground surface but not from the line itself.

I want to get it done in ArcGIS.

I am coding in VB.NET language

I am adding a rough 2-D representation to understand it better


texts:
How can I use a DBF layer in an R script in QGIS?

For example, I can input a vector by using ##Layer=vector, and a CSV table with ##Table=table, but it will fail if the table is based on a DBF file.

I can export from DBF to CSV, but then I will loose the joins I have defined.


texts: find polygons within a radius of one point (of a table of points)

texts:
I'd like to perform a select by location analysis in PostGIS that first selects the geometry of ONE point out of a table (of many points). Creates then a circumference around it and checks which geometries (polygons) of a second table lie within that circumference.

I came up with this but am somehow stuck..

SELECT t2.name, t2_area
FROM ( SELECT geom
       FROM views.sites
       WHERE site_id LIKE 'AB_10003' ) t1, spatial.land_use t2
WHERE ST_DWithin(t1.geom, t2.geom, 1000) -- radius of 1000m

the SRID of my geometries is 31467. Any help?


texts:
I've got a typical two-dimensional array in Java with values which represent a grid. Rows, columns, coordinate system (latlon), origin and cell size are known. I am able to export this data as ASCII Grid but I need to create a compressed GeoTiff. And I would like to stay in the Java world, so doing it with gdal_translate is straight forward but not the solution I am looking for.

I found GeoTools GeoTiffWriterTest.java but am not able to extract the main part I need.


texts:
I'm trying to set the symbology of a raster output in ModelBuilder of ArcGIS 10.1.

I've followed the guidelines by creating the .lyr file and the adding it to the Layer Symbology tab of the Output properties.

The model works fine but for some reason the output will not take in the style from the pre-defined .lyr file and the output still have default symbology when I run it.

Is there anything else missing here?


texts: Changes in "processing" in qgis 2.18.10?

texts:
I've been runing following code:

processing.runalg("grass:v.voronoi", self.input_file_path + file_name + ".shp", False, True,
                      extent, -1, 0, 3, file_name_poly)

or an example in plain text:

processing.runalg("grass:v.voronoi", 'C://Users/Matrix//Downloads//MsaviData//AgData//Prescriptions//MsaviData.shp', False, True, '13.539231,13.547340,55.391727,55.398065',  -1, 0, 3,'C://Users//Matrix//Downloads//MsaviData//AgData//Prescriptions//MsaviData_poly.shp')

It works (I get the output) when I run it in 2.18.0 (32bits) hence when I run it on in my 2.18.10 (64bits) it does not work? (For both test I run them in an osgeo4w installation.)

When I write "does not work", I mean it does not produce any output, if I run it in the python console it returns an output object (when I run it in 2.18.0 I got the same output message and I got the files, in 2.18.10 just the output message ). There is no error message etc.

Any ideas?
Best regards
Axel


texts: Stacking and masking Landsat at the same time with rasterio

texts:
Based on answers from previous questions 1,2 and the Masking a raster using a shapefile guide from rasterio cookbook, I wrote the code below.

The code runs correctly. I would like to improve it.

1- Notice that I'm using rasterio.mask.mask(src0, features,crop=True) twice. is there a way to keep only one?

2- Before the last line of code I'm slicing out_image[0,:,:] array in order to get a 2D array. Can avoid this operation?

import rasterio
import rasterio.mask
import fiona

# rasterio.__version__ =='1.0b1'

#open a shapefile, read geometries
with fiona.open(shapefile, "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]

# Read metadata of first file
with rasterio.open(list_to_stack[0]) as src0:
    meta = src0.meta

    # extract new mask raster info
    out_image, out_transform = rasterio.mask.mask(src0, features,crop=True)

# Update meta due to stacking
meta.update({"count": len(list_to_stack), #set the number of stack layers
             "driver": "GTiff"}) #set the input is ENVI format

# Update meta with masking info
meta.update({"height": out_image.shape[1],
             "width": out_image.shape[2],
             "transform": out_transform})

#Read each layer and write it to stack
with rasterio.open(output_stack, 'w', **meta) as dst:
    for id, layer in enumerate(list_to_stack):
        with rasterio.open(layer) as src1:
            out_image, _ = rasterio.mask.mask(src1, features,crop=True)

            dst.write_band(id + 1,out_image[0,:,:].astype(rasterio.int16))
            # add source dataset description
            dst.set_band_description(id+1 , src1.descriptions[0])


texts: Lambert_Conformal_Conic Geographic coordinate system

texts:
A shape file I am working on uses the Lambert_Conformal_Conic projection.  I would like to get a reference to this system via the IGeographicCoordinate interface.

To get at the standard lat long projection, I do this

IGeographicCoordinate gcs = srFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);

So I looked for lambert in the enum type of esriSRGeoCSType but couldn't find anything close. Does anyone know what I'm missing?

I'm using ArcEngine 10 in VS2010 with C#


texts: Exporting separate images from one layer with categorised data

texts:
this may be a poor first question but it's been niggling me most of the afternoon!

I'm trying to export a series of images from QGIS 3 based on the categorisation of the data in one layer but only show one category at a time. I've done similar before in MapInfo but that meant copy/paste and then altering the query in a lot of text files in order to isolate each category before exporting it out by hand. Is there a way of doing this within QGIS in order to automate it?

I keep thinking I'm overlooking something basic here in order to select and single out one category at a time and export the image.


texts: QGIS combine attribute table of different layers

texts: Is there a way to create a single table that contains multiple layers with all of their attributes? I have multiple layers that have the same feature options (but with different values). I basically want to add rows onto the table that come from different layers. Can this be done? I have it so my individual layer attributes automatically update as changes are made to the layer, and I would like automatic updates to show up in this new combined table as well. I have example attribute tables of two of my layers below.


texts: Set "MultiPatch" as feature class filter in arcpy scripts (ArcMap10)

texts:
As mentioned in ArcMap10.X documentation (http://desktop.arcgis.com/en/arcmap/10.3/analyze/creating-tools/programming-a-toolvalidator-class.htm), it is possible to set filters in tool parameters of python tool scripts.

So for example, for a parameter type "Feature Class" or "Feature Layer", the following geometry filters are selectable from a check list:

  A list of allowable feature class types, specified with the values
  "Point", "Multipoint", "Polyline", "Polygon", "MultiPatch", "Sphere",
  "Annotation", "Dimension". More than one value can be supplied to the
  filter.

So, when selecting a layer or feature class, only the files with selected geometries are shown.

However, when I try to set filter, the option "MultiPatch" doesn't appear in the checklist (see image below):

According to documentation, "Multipatch" and "Sphere" options should appear in the checklist, but they don't. Does anyone have an idea of why this is happening?. Another solution would be programming a ToolValidator class to set the "MultiPatch" option, using guidelines in the link above, but I don't know exactly how to proceed in this case.


texts:
How can I divide a polyline at a specific node with QGIS (version 2.12.3-Lyon)?

I was looking for something similar to MapInfo tool "Polyline split at node" which is this simple:

Make your layer editable
Select the polyline to be splitted
Click on "Reshape" (which shows all nodes)
Select the node where you want to split the polyline
From the menu "Objets" select "Polyline Split at Node"


texts: Why isn't this Python calculation working correctly?

texts:
I want to label the source files for a unified shapefile, and there are three respective columns for each. If there is no value in the column, there is no data from that source; if there is a value in multiple columns, it is from multiple sources (and therein lies the problem). I have written the following for use in the Codeblock in field calculator:

def srcfind (grcode,grcode2,lur):
 if lur == 'Parks and Recreational' and grcode == 0 and grcode2 == 0:
  return 'dmti'
 elif lur == 'Parks and Recreational' and grcode > 0 and grcode2 > 0:
  if grcode <> 0 and grcode2 <> 0:
   return 'dmti, land cover, clutter'
  elif grcode == 0 and grcode2 > 0:
   return 'dmti and land cover'
  elif grcode > 0 and grcode2 == 0:
   return 'dmti and clutter'
 elif lur <> 'Parks and Recreational' and grcode > 0 and grcode2 == 0:
  return 'clutter'
 elif lur <> 'Parks and Recreational' and grcode == 0 and grcode2 > 0:
  return 'land cover'
 elif lur <> 'Parks and Recreational' and grcode > 0 and grcode2 > 0:
  return 'clutter and land cover'

The first 'if' works correctly; even the 'dmti, land cover, clutter' category works correctly; and the last three elifs work correctly, it's just the middle part to determine 'dmti and land cover' and 'dmti and clutter' which don't work, resulting in the field being labelled 'land cover' or 'clutter' (depending on which column has a value).

Essentially, how do I fix this?

The image below is the attribute table with the variables labelled. You can see that it doesn't work for the conditions I mentioned.


texts: Writing if statement into ArcGIS Online Operations Dashboard / list options to configure style according to attribute value?

texts:
How do I write if statement into GIS operation dashboard / list options to configure style according to attribute value ?

i need to give a different background color to the list results according to attribute value.


texts: Adjacency Matrix from linestring using GeoPandas

texts:
I got my adjacency matrix wrong way round.
I got the intersection of edges. I need to have information per node instead.
as per https://math.stackexchange.com/questions/1890620/finding-path-lengths-by-the-power-of-adjacency-matrix-of-an-undirected-graph

I want to use adjacency matrix to the power of x to assess a number of paths.

from shapely.geometry import *
import geopandas as gpd
import pandas as pd
import os
import descartes # package allows to plot
import matplotlib.pyplot as plt
import fiona

Crewe = road_network
df = Crewe['geometry']

#Crewe.plot(alpha=1)

#Adjecency Matrix that takes value of 1 when there is intersection and 0.0001 if there is not.

line_gdf = Crewe.iloc[0:2500]
ma = line_gdf.geometry.apply(lambda g: line_gdf.intersects(g))
lll = ma < 1  # Where values are low
ma[lll] = 0  # All low values set to x
ma2 = ma
ma2

#Generating Adjecency Matrix to the power of x
from numpy.linalg import matrix_power
ma3 = matrix_power(ma2, 2) #instead of 3 we can take any number
ma3

My matrix looks like that now:
I want 0 to represent a node. Now I have 0 as an edge and I know it interacts with edge 4.

0   1   2   3   4   5   6   7   8   9   ... 2340    2341    2342    2343    2344    2345    2346    2347    2348    2349
0   1.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1   0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2   0.0 0.0 1.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3   0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4   1.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2345    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0
2346    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0
2347    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0
2348    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0
2349    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0


texts: assigning random points to unique values

texts: I have generated random points in qgis and would like to assign the value of their geometry to a unique column in another table. I would like to do this in qgis or postgres or spatialite. Table A contains the random points and geometry. Table B contains non spatial data and an atrribute ogc_fid which is unique. I would love to assign for each unique ogc_fid the geom of one particular feature in Table A. Any ideas of how to do this. I am comfotable with sql