texts: Trouble importing NumPy in GRASS. Architecture

texts:
I can't get GRASS to work on my system, and after doing some digging I think it is related to me using 64-bit python. I'm not really sure. I get this message...

        This module requires the Numeric/numarray or NumPy module,
        which could not be imported.  It probably is not installed
        (it's not part of the standard Python distribution). See the
        Numeric Python site (http://numpy.scipy.org) for information on
        downloading source or binaries.

when I try to import numpy from the GRASS shell I get:

ImportError: dlopen(/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-
packages/numpy/core/multiarray.so, 2): no suitable image found.  Did find:
    /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-
packages/numpy/core/multiarray.so: mach-o, but wrong architecture

I am using GRASS 7.0.0. I installed all of the required frameworks without a problem. I currently have python27-apple selected, but I have other versions on python from macports. Numpy imports fine into any version of python I select, but not into GRASS. Can anybody help?


texts: Importing R-generated .gri/.grd raster files in GRASS?

texts:
I received data in .grd + .gri format, where the first one is a small text file and the second one (having the same name) is a binary file.

I am unable to open them in GRASS/QGIS. Any hint ? Missing library ?
I tried to use Google but unfortunately "r" word is too small, and it ends to refer always to the r.* commands in GRASS...

This is the content of the .grd file:

[general]
creator=R package 'raster'
created= 2013-12-14 15:44:30
[georeference]
nrows= 134
ncols= 143
xmin= 56000
ymin= 1613000
xmax= 1200000
ymax= 2685000
projection= NA
[data]
datatype= FLT8S
byteorder= little
nbands= 5
bandorder= BIL
minvalue= 0:0:0:0:0
maxvalue= 0.945:1:1:1:1
nodatavalue= -1.7e+308
[legend]
legendtype=
values=
color=
[description]
layername= Value:Twin:Psum:pH:CN
history=


texts: Getting bar chart with negative values to display in ArcGIS for Desktop?

texts:
Using ArcMap 10.4 I'm trying to display bar charts to represent the changes that occurred at several locations. At one point, the change was always negative, and so are the three values (-3; -1; -4). For whatever reason, the chart at this point won't display.

I think it has to do with the values all being negative. There is no issue when negative values are displayed with a positive value.

Any suggestions on how to get a bar chart with negative values to display?


texts: Dissolve Leaflet VectorGrid polygon at tile boundary by ID

texts:
I am using Leaflet and the Leaflet.VectorGrid plugin to display vector tiles on a map. My vector tiles contain Polygons that render fine, except that they are also rendering the split at the tile boundary (see image below).

I have added the getFeatureID function to the creation options, so that the VectorGrid plugin knows which parts of a polygon belong together.

Is it possible to not draw the polygon boundaries at the edges of the tiles, so that it looks like one large polygon?

I did not find out how to do it (other than using a fill sytle but no outline). If it is not possible with this plugin, is there an alternative?


texts: Creating Shapely LineString from two Points

texts:
If have two points, from which I want to create a straight LineString object:
from shapely.geometry import Point, LineString
A = Point(0,0)
B = Point(1,1)

The Shapely manual for LineString states:

A sequence of Point instances is not a valid constructor parameter. A LineString is described by points, but is not composed of Point instances.

So if I have two points A and B, is there a shorter/better/easier way of creating a line AB than my current "best" guess...
AB = LineString(tuple(A.coords) + tuple(B.coords))

... which looks rather complicated. Is there an easier way?

With Shapely 1.3.2, the above statement from the manual is no longer correct. So from now on,
AB = LineString([A, B])

works!


texts: I have some prj files (Only prj) but don't know how to open them, I tried ArcMap but have no clue. I understand prj is a part of a shapefile, but I don't have any other files than these.


texts: Convert an ASCII grid file to GeoTIFF using Python?

texts:
I have a ASCII grid raster format file. For example:

ncols 480
nrows 450
xllcorner 378923
yllcorner 4072345
cellsize 30
nodata_value -32768
43 2 45 7 3 56 2 5 23 65 34 6 32 54 57 34 2 2 54 6
35 45 65 34 2 6 78 4 2 6 89 3 2 7 45 23 5 8 4 1 62 ...

How can I convert it to TIFF or any other raster using Python?


texts:
Is there a method in Arc Gis for mapping storm water runoff in an urban environment assuming that the sewers are at full capacity. I just want to map roughly where the water will flow over the land over paved surfaces

I have a DEM, a land use map and have already calculated the flow direction and streams network. I wish to find out what areas will flood in the study area during an extreme rainfall event.


texts:
I'm trying to interpolate the depths of ground water aquifers with a few known data points from the field study. The objective here is to map the contours of that particular catchment's depth of bedrock or aquifers so that we can strongly assume the ground water depth values for the ones we don't have a depth value?
 Our study area is on the peri urban small towns of Bangalore, India


texts: Convert KML to GeoJSON with ogr2ogr and remove KML special attributes

texts:
I'm trying to convert KML to GeoJSON (or some other formats) using ogr2ogr  (look at this question). It's all working fine, except of one problem:

I'm always receiving these  metadata attributes with any operation (using org2ogr or ogrinfo or Python script)

  Name, description, timestamp, begin, end, altitudeMode, tessellate,
  extrude, visibility, drawnOrder,icon

But my Schema doesn't have these fields:

<Schema name="distritos_mad" id="distritos_mad">
    <SimpleField name="ID" type="float"></SimpleField>
    <SimpleField name="CD_GEOCODS" type="string"></SimpleField>
    <SimpleField name="NM_SUBDIST" type="string"></SimpleField>
    <SimpleField name="Nomes" type="string"></SimpleField>
    <SimpleField name="AREA (HA)" type="float"></SimpleField>
    <SimpleField name="Pop (2010)" type="float"></SimpleField>
    <SimpleField name="Ddem p/km" type="float"></SimpleField>
</Schema>

Looking at LIBKML Drive documentation I realize that this fields are "special fields" and for an unknown reason this fields are included on feature properties.

How can I prevent OGR from writing these feature attributes?

I'm runing this command:

ogr2ogr -f "GeoJSON" C:\out.json C:\in.kml -t_srs EPSG:4326 -dim 2 --config OGR_FORCE_ASCII NO

  GDAL 1.11.3, released 2015/09/16


texts: Changing from layout view to data view using ArcPy?

texts:
Is there a way to change from layout view to data view using Python and ArcMap 10?

I would prefer to do it using ArcPy.


texts: How to Add (+) up integers in a field and print a total count in Python 2.7

texts:
I'm a complete beginner to python. I have what I hope to be a simple question. I want to run a bit of script in IDLE so that it adds up all figures in a table under attribute name "Tot_Aff". When I run the below it just gives me each Value. Any help would be great

fc = final_SOA
field = "Tot_Aff"
cursor = arcpy.SearchCursor(fc)
for row in cursor:
    print row.getValue(field)

I want to add up all the integers


texts: Can't add POINT NOT NULL column to MySQL table

texts:
I recently upgraded to MySQL 5.7.9 & am using MyISAM for my table engine.

Adding a new POINT column like so:

ALTER TABLE <table> ADD centroid POINT NOT NULL;

gives this error:

  ERROR 1138 (22004): Invalid use of NULL value

I can create the column without NOT NULL but then I am unable to create a spatial index on the column. This worked for me in MySQL 5.6, and I don't see anything relevant in the docs.


texts:
I am trying to set up 3d city database and I have some problems. I hope you can help me. I receive an error as you can see attaching file. What is wrong?

The error message is:

'"\psql"' is not recognized as an internal or external command, operable program or batch file.

and results from this command:

C:\long_path\PostgreSQL\SQLScripts>"\psql" -d "3dcitydb PGBIN=C:\Program files\PostgreSQL\9.5\bin" -f "CREATE_DB.sql"


texts: Referencing a Numeric Value in ArcGIS ModelBuilder?

texts:
How do you reference a numeric User input number in ModelBuilder?  I built a ModelBuilder tool that creates a traffic map based on a users input.  One question is between what two speeds is the traffic going?  So I need to select everything above one speed and below another......

In Select by attribute I tried Speed>%MinSpeed% but it did not work.

Both the user input variable and field type are 'double'.

The one possibility that just occurred to me is the numbers I have been using as user input numbers are whole numbers however the real data has it down to the decimal place. Would that effect anything? So user input is 20 and real data could be 20.05.


texts: Tile Package made with ArcPy are not tiling all layers or levels of detail

texts:
I'm having a weird issue with CreateMapTilePackage_management in ArcPy. The tile packages that are created by this tool do not have all layers and levels of detail. I'm not sure what is causing this. Even the tiles look off when put into arcmap.

Has anyone ever seen this?

Here is the code example:

 df = arcpy.mapping.ListDataFrames(map_document, "Road_Map")
 for lyr in arcpy.mapping.ListLayers(map_document, "", df):
    if lyr.name == "Conductor":
       print "\n\n**** FOUND THE CORRECT MAP LAYER ****\n\n"
       time.sleep(3)
       print "\n\n**** SETTING UP ARCGIS WORK ENVIRONMENT ***\n\n"
       arcpy.env.workspace = loc_of_featureclasses
       for fc in arcpy.ListFeatureClasses():
            print "\n\n**** CHANGING LAYER DATA SOURCE ****\n\n"
            lyr.replaceDataSource(loc_of_featureclasses, "FILEGDB_WORKSPACE", fc, True)
            map_document.save()
            print ("\n\n**** ZOOMING TO LAYER EXTENT ****\n\n")
            dataFrame = arcpy.mapping.ListDataFrames(map_document)[0] # the first data frame
            MapLayers = arcpy.mapping.ListLayers(map_document,"Conductor",dataFrame)
            Layer = MapLayers[0]
            Extent = Layer.getExtent(True) # visible extent of layer
            dataFrame.extent = Extent
            arcpy.RefreshActiveView() # redraw the map
            time.sleep(3)
            print "\n\n**** SETTING UP ARCGIS WORK ENVIRONMENT ****\n\n"
            arcpy.env.overwriteOutput = allow_tile_to_overwrite
            arcpy.env.workspace = tile_workspace
            time.sleep(5)
            print "\n\n**** CREATING TILE PACKAGE FOR MAP " + lyr.datasetName + " ****\n\n"
            arcpy.CreateMapTilePackage_management(mxd_file_name, "ONLINE", lyr.datasetName + ".tpk", "PNG8", "10")
            print "\n\n**** TILE PACKAGE FOR MAP " + lyr.datasetName + " FINISHED ****\n\n"
            time.sleep(6)
            print "\n\n**** MOVING ON TO THE NEXT MAP ****\n\n"


texts: Digitize by User Input by Distance and Length in QGIS

texts:
I was wondering if anyone knows a way to digitize a line or polygon by user entered distance and length.  This is an option in ArcMap, but I would like to be able to do this in QGIS.

Thanks,
Ryan


texts: How to display shape files stored in a database into a web page

texts: I'm pretty new to geodjango and python in general. But I would like to display my shapefiles stored in my database into a web pages and present them in layers and retrieve certain information if a user want to.


texts: How to create vector polygons that show areas of overlap (2, 3, 4 times, etc.)?

texts:
I have polygon shapefile containing several thousand features.  The attribute table contains, among other fields, year (1943-2019).

https://imgur.com/a/C2JTIjX

These polygons (wildfires) often overlap.  I would like to create a shapefile that inputs this layer and outputs a polygon shapefile of all areas that overlap twice, three times, four times, etc.

I would be especially happy if I can have the years of the polygons that overlapped included in the table of these proccessed rasters.

I can also parse the wildfires polygon into separate shapefiles based on the year of the fire.

I will prefer to use ArcGIS Pro 2.4.  I can also use ArcMap 10.7 or QGIS 3.8.3 if your proposed solution works for a specific program.


texts: How to export dynamic track (video) from QGIS?

texts:
Is it possible to export a file video with QGIS with a dynamic track from point A to point B?

I need something like a presentation of a race where a dot is moving along the track from the start to the finish

I use QGIS 2.6.1


texts: Setting the extents or bounds of "map view" of OpenLayers does not open the webpage at given Zoom Level

texts:
When I load the map without any extent in my view, then the HTML page opens properly,

      var myview = new ol.View({
            projection: 'EPSG:4326',
            center: [77, 24],
            zoom: 4.5,
            maxZoom: 8,
            minZoom: 3,
            smoothResolutionConstraint:true,
            //extent: [60, 6, 100, 40],
      });

but once I set the bounds in my view, it doesn't follow the given zoom level and opens differently

      var myview = new ol.View({
            projection: 'EPSG:4326',
            center: [77, 24],
            zoom: 4.5,
            maxZoom: 8,
            minZoom: 3,
            smoothResolutionConstraint:true,
            extent: [60, 6, 100, 40],
      });

How should I change the JS so that, it opens at the given zoom level when I give the bound extents.


texts: Display extent indicators in the legend?

texts: Can you link the data frame extent indicator box with your legend so that the legend displays a symbol and description representing the box?


texts:
I recently posted my first question on StackOverflow, but I think the question was quite extensive and full of potential errors. So I'd begin again, step-by-step, with a much shorter and simplier question :

The purpose is to structure my models.py in my django website to store datas in a Mongo database as GeoJson objects.

Do you think this code is correct, I'm not very sure about the syntax ?

Thanks a lot !

models.py :

# -*- coding: utf-8 -*-

from mongoengine import *
connect(‘mongodb_jsons’)

import datetime

class GeoJson(Document):
 # from a geojson object
 # save it in Mongodb with a geojson structure

    Date_created = DateTimeField(default=datetime.datetime.now)
    Latitude = FloatField()
    Longitude = FloatField()
    Content1 = CharField()
    Content2 = CharField()

    meta = {'db_alias': 'mongodb_jsons', # save in DB ‘mongodb_jsons’
            'indexes': [                 # the geojson structure
                {'type' : 'Feature', {
                   'geometry':{
                        'type' : 'Point',
                        'coordinates' : ['Longitude','Latitude']},
                        },
                   'properties':{
                       'content1' : 'Content1',
                       'content2' : 'Content2',
                       'date_creation' : 'Date_created'
                        }
                    }
                }]
            }


texts: Displaying label outside canvas/map extend in QGIS

texts:
The first picture shows my map with cadastre (black lines), points with label (orange squares/ red numbers) and a printframe (blue rectangle). The second picture shows what I want in print layout.  It consists of two rectangles, one with the printframe from the map and one with the labels (without map and for now only with manually places text boxes). I want to label the shown points in the map in a second map above (or any other way). The labels need to be placed on the same vertical and horizontal line. I tried to display the labels on x,y (x-fixed for horizontal line), but couldn't get them displayed because items are not in the map/canvas and won't be shown. Furthermore I had some difficulties because of the rotation of the map/printframe.
One solution might be to shown the labels even when the objects are not in the canvas/map extend. Is that somehow possible?
The only other solution I can think of is to get the x,y-values of every point (>500) within every printframe (>100) I have, add a individually x,y-value to every point based on the height of the frame for the labeling and have two maps in the print layout on the top of each other.
As I don't want to do it the other way, it would be nice if someone can give me a hint.


texts: ERROR Cutline feature without a geometry

texts:
I am using QGIS 2.18.24. I'm trying to clip raster (asc) by mask (shp).
I'm working in EPSG:4326.

I'm checking the geometry of the polygon file, no errors were found.
Geometry is valid.
Edit mode is disabled.

I'm always getting the following error message that I haven't been able to sort out:
Algorithm Clip raster by mask layer starting...
GDAL command:
gdalwarp -ot Float32 -q -of GTiff -cutline C:\Users\Urbana\Desktop\tesi\Europa_shp\EUROPA_OK.shp -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -wo OPTIMIZE_SIZE=TRUE "C:\\Users\\Urbana\\Desktop\\tesi\\variabili current\\Nuova cartella\\bio01.asc" C:\Users\Urbana\AppData\Local\Temp\processingf8e7cc09212f4e19a00d49110f0c0b1c\a879020e1c284cfaabb238991ac727b1\OUTPUT.tif
GDAL command output:
ERROR 1: Cutline feature without a geometry.
Converting outputs
Caricamento layer

The following layers were not correctly generated.
Ritagliato (da maschera)
You can check the log messages to find more information about the execution of the algorithm


texts: Reprojecting multiple rasters in ArcGIS Desktop?

texts:
I have classified layers of Landsat 5 and 7 of several years in Geographic coordinate system WGS84.

To reproject this to UTM, is there a simpler way than using Project Raster tool in ArcToolbox?


texts: Testing Validity of Web Service Referenced in Layer File?

texts:
My organization creates a number of layer files for our ArcMap users to access data from web services. Is there any way to use arcpy.Describe to get to the URL of the service referenced by the layer? My testing so far has only yielded paths back to the layer file itself, or a temp folder:

>>> d=arcpy.Describe(r'e:\geodata\ortho_imagery\naip_ca_historical.lyr')
>>> d.catalogPath
u'e:\\Geodata\\ortho_imagery\\naip_ca_historical.lyr'
>>> l=d.layer
>>> l.catalogPath
u'C:\\Users\\bjorn\\AppData\\Local\\Temp\\3\\NAIP_Historical/California_Historical'

I'd also be open to a .NET solution if necessary.

The ultimate goal is to develop a tool that will scan our shared drive for .lyr files, and if they refer to a web service, test its validity.


texts:
i work with RStudio and Leaflet. I like to add inputboxes to type coordinates into them. Then the code should set a marker on this point and zoom on it.
How can i do this?
I tried to manipulate the html code with: but it does not run :( . I think there is a missing specification on my current map.
Please help :)

      X:
      <input type="text" id="x">
      Y:
      <input type="text" id="y">
      <button type="button" onclick="goThere()">Go</button>

var secondProjection = "+proj=lcc +lat_1=29.7 +lat_0=29.7 +lon_0=-5.4 +k_0=0.9996155960000001 +x_0=500000 +y_0=300000 +a=6378249.2 +b=6356515 +towgs84=31,146,47,0,0,0,0 +units=m +no_defs";
  var firstProjection ="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "


texts: Using global parameters in nested models of ModelBuilder?

texts:
I have a model that consists of two nested models and various parameters. There are two parameters that are common to all of the models (including the overall “grand” model). The common parameters are in-line variables and serve to assign unique file names to the models’ output.

Is there a way to set the values of these parameters at a single time and then have these values propagate through the nested models?


texts: ArcGIS REST Services Directory url without the port 6080

texts:
I am wondering if I could access the GIS server URL http://gisserver.domain.com:6080/arcgis without the port(6080) such as http://gisserver.domain.com/arcgis. I was told that Using a reverse proxy server with ArcGIS Server. I am wondering if anyone could provide your experience how to hide the port.

The background of my server

ArcGIS server 10.3.1 (No Web Adaptor installed)

I have Windows 2012 Server x64 with SQL Server 2012.

I'm also using IIS 8 as a web server on the same machine.

It is the security issue. In my case, I have enabled the windows authentication for the website "gisserver.domain.com";. However, the window authentication does not work when I type the "gisserver.domain.com:6080";. I have set the inbound rule in the firewall setting inside the server, which opens the port 6080 and 6443 to the public. I was always told to use the Web Adaptor to address this issue. I am wondering if there are other methods to address this issue in the IIS or firewall setting.


texts:
So I have the following command line that runs in my php code
First:
I rename the layer_name with the id of the user who uploads his shp file,
and save it with the kml extension
ogr2ogr -f KML  11141601905540.kml  141601905540.shp   -nln 141601905540

Second:
a).-I connect to my postgis database with postgress and save the kml information.
b).-Change the coordinate reference system from EPSG:4326 to EPSG:3857:
ogr2ogr  -f PostgreSQL -s_srs EPSG:4326 -t_srs EPSG:3857    PG:"host=127.0.0.1 dbname=postgis_24_sample port=5432 user=postgres password=123456"   C:\xampp\htdocs\siap-padron\public/uploads/AAGL770425MGRYNL05/zip/shp/11141601905540.kml  -nln  layers_fisico  -dialect sqlite  -sql "SELECT  '11141601905540' as id_unico, 'LILIA' as name,'1' as category_id, '3' as food_industry, 'AAGL770425MGRYNL05' as curp,'03' as cve_edo,'001' as cvemun,'69800' as zip_code,'0008' as location,  GEOMETRY AS wkb_geometry  FROM \"11141601905540\""  2>&1

The code works great, until one of the users couldn't upload his shp file
so I put the content with file A (IT IS LOADED WELL) and a file B (IT IS NOT LOADED)
A)
<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document id="root_doc">
<Schema name="151601905804" id="151601905804">
<SimpleField name="Id" type="int"></SimpleField>
<SimpleField name="RazónS" type="string"></SimpleField>
<SimpleField name="Cul_Esp" type="string"></SimpleField>
<SimpleField name="RFC" type="string"></SimpleField>
</Schema>
<Folder><name>151601905804</name>
 <Placemark>
   <ExtendedData><SchemaData schemaUrl="#151601905804">
    <SimpleData name="Id">4</SimpleData>
    <SimpleData name="RazónS">Agropecuaria del Sur</SimpleData>
    <SimpleData name="Cul_Esp">Bovino-carne</SimpleData>
    <SimpleData name="RFC">HJU890108K56</SimpleData>
   </SchemaData></ExtendedData>
  <Point><coordinates>-98.5500990493249,18.6732116243429</coordinates></Point>
</Placemark>
</Folder>
</Document></kml>

B)
<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document id="root_doc">
<Schema name="151601907584" id="151601907584">
<SimpleField name="RFC" type="string"></SimpleField>
<SimpleField name="RAZONS" type="string"></SimpleField>
</Schema>
<Folder><name>151601907584</name>
<Placemark>
<Style><LineStyle><color>ff0000ff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>
<ExtendedData><SchemaData schemaUrl="#151601907584">
    <SimpleData name="RFC">HGH851CDCD</SimpleData>
    <SimpleData name="RAZONS">Sin nombre</SimpleData>
</SchemaData></ExtendedData>
  <Polygon><outerBoundaryIs><LinearRing><coordinates>0.0,821216.592200002 0.0,821290.949 0.0,821466.091399997 0.0,821395.670100002 0.0,821216.592200002</coordinates></LinearRing>
</outerBoundaryIs></Polygon>
</Placemark>
</Folder>
</Document></kml>

file B generates the following output
ERROR 1: Failed to reproject feature 0 (geometry probably out of source or destination SRS).
ERROR 1: Terminating translation prematurely after failed

I have seen that there are many projections, so how could I solve this?
Any comment or article to investigate will be helpful


texts: Difference of georeferenced and ungeoreferenced map in AutoCAD?

texts: Can somebody please explain the difference of an ungeoreferenced scanned map and a georeferenced one in AutoCAD?


texts: Flex Map component, update extent on resize

texts: I have application with Flex Map component. Application has several states, and map component's size is changed depending on states. When state is changed, I need to change extent. But for some reasons, map component is not recalculated by that time. I tried to change extent on resize event. I see in debug mode, that map width is changed, but when I try to set new extent, it "thinks" that map width is not changed yet. So, every time I set new Extent, it's set for the previous map state width. I tried callLater, but all the same, map component has not calculated something inside itself and sets extent for the previous map width. Thanks in advice


texts:
I am trying to cache my raster tiles from Tileserver-GL. This is my actual link of tiles
http://xxx.xx.xx.xx:8080/styles/map/{z}/{x}/{y}.png
(Each tile link will look like http://xxx.xx.xx.xx:8080/styles/map/1/1/1.png)
I have to change the default.vlc file to
backend default {
    .host = "127.0.0.1";
    .port = "8080";
}

Now I am getting the tiles from the following link
http://xxx.xx.xx.xx:6081/styles/map/{z}/{x}/{y}.png
(Each tile link will look like
http://xxx.xx.xx.xx:6081/styles/map/1/1/1.png)
When I am doing my first map load varnish is calling the tiles from the Tileserver-GL server. But the problem, it is still calling the tiles from the Tileserver-GL server. That means there were no cache png files stored.
Here is an example of how OpenStreetMap is caching there tiles. They have 3 layers but in my case I have only one layer.
https://wiki.openstreetmap.org/wiki/Varnish_setup


texts: Iterate through multiple Raster files and compute statistics using Zonal Statistics

texts:
I have a large number of raster images. What I'd like to do is to compute statistics on different sections of these raster files using a polygon. I think everyone is suggesting Zonal statistics. However, I'd like to iterate through all the raster datasets and then store their statistics in a single table. Zonal statistics does this for one raster and then overwrite it next time. I used the %Name% inline variable as for my table name, but that just creates lots of tables.

How can I iterate through all this raster datasets and create only one table?


texts: Defining a line from geographic coordinate points.

texts:
On a flat surface we can define any point on line which passes from point1 (x1,y1) to point2 (x2,y2) via simple trigonometry.

How can we do that with latitude and longitude of 2 points?

Imagine we have geographic points (54.977614, -103.710937) and (55.973798, -107.929687) on a globe. Using what mathematics can I define any point on the line which passes through these points? Or where can I read about this?


texts: How to make GRASS and SAGA tools capable to use in QGIS？

texts:
I installed QGIS, GRASS and GDAL by OSGeo4w package, and as I finish the installation, I want to use GRASS tool, I face this problem：

But I already installed GRASS and SAGA, how to make those tools in QGIS to be used？


texts: Getting lat/lng from wkid, latestwkid and x y coordinates

texts:
I am receiving some JSON feed with this data
"geometryType": "esriGeometryPoint",
"spatialReference": {
    "wkid": 102100,
    "latestWkid": 3857
  },

and set of coordinates like
"geometry": {
   "x": -11696523.780400001,
   "y": 4804891.0001000017
}

If there a way to convert this x and y to lat/lng which I can use in Google Maps?
I have read some references but there either I can`t see the answer, or it's goes very deep into GIS terminology which I just can't afford to learn at the moment.


texts: ArcGIS 10.1 Crashes when updating address locator

texts:
I keep getting a "ArcGIS has encountered a serious error and needs to close" crash dialogue box whenever I try to update an address locator imported from one server to another (updating the location of the reference database).

I've tried closing ArcMap and just using ArcCatalogue, still doesn't work.

I'm on a 3Ghz/4G RAM relatively new (<1 year) PC at work.

Can anyone help me with a work-around for this?


texts: Modifying a map object parameter for ArcGIS JavaScript API printTask?

texts:
I want to change the extent and operational layers of a map for several printTasks. The goal is to take a group of polygons, and create a map of each one separately, changing the basemap and some operational layers on each map. The resulting maps will end up in a report. What is the best approach to do this? A couple thoughts I had:

Change the map object to the desired operational layers and extent before printing, then revert the map object back to the original state after printing. (seems like a lot of extra steps this way)
Use a modified Web_Map_as_JSON parameter for each map. I could have this is a pre-configured string that would be easily modified. If this is a reasonable solution, is it possible to add Web_Map_as_JSON` as a parameter to the printtask like this:

    //pre-configured JSON map string
    var web_map_as_json = this.web_map_as_json

    //modify JSON here (extent of each polygon, etc..)

    var params = {
        Web_Map_as_JSON: web_map_as_json,
        Format: 'png32',
        Layout_Template: 'MAP_ONLY'
    };

    printTask.execute(params, lang.hitch(this, 'onParcelPrintComplete'), lang.hitch(this, 'onParcelPrintError'));


texts: Transforming / stretching map or it's part in particular direction

texts:
I am trying to model ancient land surfaces we can 'rollback' geological folding and faulting. This is called palinspastic modelling or mapping.

So I have an area that due to folding and faulting was reduced by 50% in a N-S direction. I am using QGIS.

I have been unable to find any preset function that will allow me to stretch my map by 150% in the NS direction.

I did look at creating a custom CRS but even within proj and WKT there does not appear to be this facility. I am sure there must be and I have missed it.

Can anybody help?


texts:
I georeference a tif image and then I clipped that image in a mosaic.

I did a index of these images but when I load them is very slow(7 seconds or more to display the first tile). How can I do to go faster?

I tried many ideas from the mapserver document but I does't work. Any idea?


texts: Trying to convert a shapefile that seems to be using UTM coordinates to WGS84 in QGIS

texts: I've been trying to convert this polygon shapefile of historic Swedish parishes from UTM coords to lat/long and nothing seems to be working. I've tried defining projection, reprojecting, saving these new projections as new shapefiles and loading them to a new project, fixing geometries, extracting vertices and converting the csv file of coords to lat/long and the closest I can get is:


texts: Unable to loop/delete with QgsVectorLayer and dataProvider

texts:
I have a csv file with the following columns:

table   geometry    field   type    typeName    len   precision

and my goal is to create a number of shp files programatically with PyQGIS.

My script stops after creating the first shp file.

Here's my code.

import sys
import re
import csv
from PyQt4.QtCore import *
from PyQt4.QtGui import *

reload(sys)
sys.setdefaultencoding('utf8')
final_list = []

with open(r"D:\Users\Ravi_narayanan\Documents\shpgen.csv", "rU") as s:
    reader_tables = csv.DictReader(s,delimiter=';')
    liste_tables = list(reader_tables)
    l_tables = [i['table'] for i in liste_tables]
    l_tables = set(l_tables)
    for table in l_tables:
        l1 = [line for line in liste_tables if line['table'] == table]
        geometry = l1[1]['geometry']
        layer_mem = QgsVectorLayer(geometry+"?crs=epsg:2154",table,"memory")
        QgsVectorFileWriter.writeAsVectorFormat(layer_mem,r"D:\setupCAPFT2.00\{}.shp".format(table),"utf-8",None,"ESRI Shapefile",False)
        layer_path = "D:\{}.shp".format(table)
        shp = QgsVectorLayer(layer_path,table,'ogr')
        pr = shp.dataProvider()
        shp.startEditing()
        for field in l1:
            if field['type'] != 'Double':
                pr.addAttributes([QgsField(field['field'],getattr(QVariant,field['type']),field['typeName'],int(field['len']))])
            else:
                pr.addAttributes([QgsField(field['field'],getattr(QVariant,field['type']),field['typeName'],int(field['len']),int(field['precision']))])
        shp.updateFields()
        shp.commitChanges()
        shp.stopEditing()
        del layer_mem,shp,pr


texts: Updating Cursor to Include Newly Created Features?

texts:
I am using ArcMap 10.1 and python to create cursor that starts at the one and only row of a polyline shapefile, some code runs to create a new feature based on the position of the exisiting feature, then adds this new feature as a second row. I then want the same code to run again on the newly created feature to then create a third...and so on.

However, the cursor stops working after the first run and the code completes, I guessing because the initial cursor doesn't take into account the newly created feature?

However I need this to run until a certain number of rows/features are created.

Is anyone able to help with this please?

Its as if I need to refresh the cursor or something, but I struggling to find a solution.

def addPolyline(cursor, array, sr):
    polyline = arcpy.Polyline(array, sr)
    cursor.insertRow((polyline,))
    array.removeAll()

import arcpy
import os
import math

infc = ...path...
shapefile = ...path...

mxd = arcpy.mapping.MapDocument("CURRENT")

df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

spatialRef = arcpy.Describe(infc).spatialReference

desc = arcpy.Describe(infc)
shapefieldname = desc.ShapeFieldName

rows = arcpy.SearchCursor(infc)

with arcpy.da.InsertCursor(shapefile, ("SHAPE@",)) as cursor:
    polylineArray = arcpy.Array()

    for row in rows:

            < SOME CODE TO CREATE NEW FEATURE >

            newStartPoint = arcpy.Point(newStartX, newStartY)
            newEndPoint = arcpy.Point(newEndX, newEndY)

            polylineArray.add(newStartPoint)
            polylineArray.add(newEndPoint)

            addPolyline(cursor, polylineArray, spatialRef)


texts:
My non-GIS coworkers are always baffled by the results of intersection queries, because they expect the result pictured below, but instead, they get every single feature being selected.

Part of the problem is that ESRI uses the term "intersect" to refer to both of these, completely different spatial relationships (see: Select by Location Intersect vs the Intersect analysis tool).

I always use work-arounds (typically a union/intersection analysis and then a series of cursors), but there must be a spatial relationship defined somewhere in ESRI's capabilities that selects features as selected in the picture below based on their relationship to the underlying, green feature. That mystery definition is apparently used in the Intersect tool, but what if we want to stop it short of cutting the features and instead copy the entire feature that otherwise would be sliced up and only partially retained?

Is there some other tool where such a selection capability might exist?

It's an "overlap", as simply and straightforwardly interpreted everyday in our language, so where is it in ArcGIS?


texts: Splitting polygon into polygons based on data point groups using ArcGIS Desktop?

texts: Is there a function in ArcGIS 10.5.1 that can split a polygon (black outline) based on the location of data points (red circles)? I need it to simply make a divide in between different groups of points, which is in this case group 1 and 2 (see image).


texts: How do I merge a point attribute table to a polygon attribute table in QGIS or GRASS?

texts:
I have a point layer containing the following data:

I want the following polygon areas to be added to the point attribute table as all the points occur inside the various features of the polygon layer:

There are thousands of points but under 50 features in the polygon layer:

PLEASE NOTE that I have attempted the following plugins/algorithms:

Merge
Join attributes by Location
Join attributes Table
spatialJoin (that uses refFunctions)

Please help!


texts: PyQGIS interpolation classes and triangulation file without interpolation raster

texts:
I am trying to interpolate lines using Interpolation plugin in QGIS. There is an option for exporting triangulation file. I don't want to create raster interpolation. In this case, I want only triangulation file. My code:

from qgis.core import QgsRasterLayer, QgsVectorLayer, edit, QgsFeatureRequest, QgsMapLayerRegistry
import qgis.analysis

resolution = 0.1
q5_edited = QgsVectorLayer("C:\\Users\\testuser\\Documents\\exp_q5_23.shp", "CONTOUR", "ogr")
q20_edited = QgsVectorLayer("C:\\Users\\testuser\\Documents\\exp_q20_23.shp", "CONTOUR", "ogr")

q5_index = q5_edited.fieldNameIndex('intpol')
q20_index = q20_edited.fieldNameIndex('intpol')

ld_q20 = qgis.analysis.QgsInterpolator.LayerData()
ld_q5 = qgis.analysis.QgsInterpolator.LayerData()

ld_q20.vectorLayer = q20_edited
ld_q20.zCoordInterpolation=False
ld_q20.InterpolationAttribute = q20_index
ld_q20.mInputType = 1 #1 = structureLines, 2 = breakLines

ld_q5.vectorLayer = q5_edited
ld_q5.zCoordInterpolation = False
ld_q5.interpolationAttribute = q5_index
ld_q5.mInputType = 1

itp = qgis.analysis.QgsTINInterpolator([ld_q20, ld_q5])
itp.setExportTriangulationToFile(True)
itp.setTriangulationFilePath('C:\\Users\\testuser\\Documents\\dualedge.shp')

rect = q20_edited.extent()
ncol = int((rect.xMaximum() - rect.xMinimum())/resolution)
nrow = int((rect.yMaximum() - rect.yMinimum())/resolution)

test = qgis.analysis.QgsGridFileWriter(itp,"C:\\Users\\testuser\\Documents\\dualedge_raster.tif",rect,ncol,nrow,resolution,resolution)
test.writeFile(True)

I set setExportTriangulationToFile to True and set the setTriangulationFilePath to path where I want to create triangulation file. But the triangulation file is created only when the last line is run test.writeFile(True). This command creates triangulation file (.shp) and raster file (.tif). I checked the code of writeFile function (https://qgis.org/api/qgsgridfilewriter_8cpp_source.html) but I can't find the place where triangulation file is created I found only the writing to a raster file. I hoped there will be some function which creates triangulation file before creating the raster file.

How the triangulation file is created and is it possible to create triangulation file without creating raster? I need it because we interpolate a very big vector data which makes raster really big and we don't need a raster, just a triangulation.


texts: Attribute Table Window loses highlighted selection after Addin is Opened

texts:
I'm new to GIS.

I have taken over a project that is an AddIn for an ArcGIS Desktop app. C# VS 2010 project. I can get the Selected Features just fine when they are selected but when I add or remove selected rows the blue highlight disappears! I know there are selected rows there as I can see the (4 out of 2000 Selected) etc at the bottom of the window. If I click show all highlighted then I see just those row which are then highlighted, and then when I click Show all Rows the highlighted row a there in blue again along with all rows, which is great.

BUT, how can I code the AddIn so that the highlights remain? Is there a Reselect Selected Features in Attribute table window method somewhere?

I declared this in my class but don't know how to use it!

public esriTableSelectionActions TableSelectionAction { get; set; }


texts: How to define voronoi polygon buffer extend in QGIS from layer extend

texts: I'm using QGIS (2.8.1) Voronoi Polygons Algorithm and I would like to extrapolate them to a polygon layer extend. In order to do this, the tool opens a window asking me to enter a number or expression (s. image) and proposing the x min/max and y min/max coordinates of all layers. What would a formular look like that covers the extend of a layer? (What do I have to enter in the field?)


texts: How to find street for point PostGIS OSM?

texts:
I need to select some information about point and street_id for this point from database. I can find minimal distance between point and street, but there are a lot of streets in result set. How can I output only one line for each point?
This is my query:

SELECT  ST_Distance(ST_LineInterpolatePoint(street_line, ST_LineLocatePoint(street_line, house_loc))
    , house_loc) As distance
    , street_name
    , street_id
    , point_name
FROM
    (SELECT line.way As street_line,
        point.way as house_loc,
        line.name as street_name,
        point.name as point_name,
        line.osm_id as street_id
    FROM planet_osm_line as line, planet_osm_point as point
    WHERE   point.amenity = 'cafe' AND
        line.name notnull AND
         ST_Contains(ST_Transform(ST_GeomFromText('POLYGON((23.692064 52.092482, 23.693534 52.093365, 23.696163 52.090194, 23.693288 52.089858, 23.692064 52.092482))',4326),900913), point.way)) as foo
WHERE ST_DWithin(street_line, house_loc, 100);

Result for first poin looks like this:

     distance     |   street_name   | street_id |      point_name
------------------+-----------------+-----------+-----------------------
 14.5456161008274 | Советская улица | 431230474 | Метромилано
 41.7027205345511 | Советская улица | 430790594 | Метромилано
 40.0583523277134 | улица Будённого | 108887692 | Метромилано

I need to get only one line with minimum distance. I tried to select min(distance) but result was the same.


texts: What is the correct way to merge adjoining raster images?

texts:
I am trying to create a raster base layer from NAIP ortho imagery. Unfortunately for me, the area of interest is very near the intersection of four image files. I have tried a few approaches to combining them into a single map layer:

Naively grouping them together in a  GeoServer Layer group
Combining them with GeoServer's ImageMosaic plugin
Merging them with gdal_merge.py

In all cases, the result has the same visual problem of ragged black lines where the images meet. See the example image, below. The black lines appear at the edges of each of the individual files, they are not merely an artifact of merging the images, so in a sense this isn't a surprising outcome -- just annoying.

What is the correct procedure to join these orthoimages into a reasonably seamless layer? Is this even a practical thing to do?

A few more details:

The imagery came from the USDA Geospatial Data Gateway, National Ag. Imagery Program Mosaic listing
The files were in MrSID format, but I've also created GeoTIFF copies for the gdal_merge.py approach
In addition to the MrSID files, the data set included shape and prj files.
The images use projection NAD83 / UTM 16N, so I've been using EPSG 26916
I'm a programmer, not a GIS expert, so my knowledge of this area is pretty limited


texts: Python GeoPandas: fill small gaps between Linestrings in a graph to get only one connected component?

texts:
I have a graph in Python GeoPandas with many connected components. And I want to create a graph with only one connected component.
Does there exist some convenience functions to fill the gaps in GeoPandas?


texts: Identify intersection points in linestrings using PostGIS

texts:
I'm trying to replicate the result of the QGIS line intersections tool in PostGIS.

The result of this tool gives me 133 intersections represented as points as depicted in the map below.

Using the same data I've tried to replicate this in PostGIS using this code

WITH select_trips AS
    (
    SELECT t.* FROM trips t
    WHERE t.id IN (25758,24730,23875,25886,24442,26723)
        )
    SELECT
        a.id,
        st_intersection(a.the_geom_webmercator, b.the_geom_webmercator) as buffer
    FROM
        select_trips a
    JOIN
    select_trips b ON st_intersects(a.the_geom_webmercator, b.the_geom_webmercator)

My problems are two fold. 1: this only returns 36 rows, 2: when I try to add the layer using the QGIS DB manager, nothing gets added.

From what I understand st_intersection should return the shared geometry of the two lines represented as a point, so why isn't this working?

I know I can run the tool in QGIS and then import to my postgresql DB, but I'd like to be able to run the whole process using PostGIS.

UPDATE:
Using this code returns some very strange results a line string that looks somewhat like the trip lines below and 1127 intersection points. (the light blue is the correct intersection points generated with QGIS line intersections tool)

SELECT
    ST_Intersection(a.the_geom, b.the_geom),
    a.id
FROM
    trips as a,
    trips as b
WHERE
    st_intersects(a.the_geom, b.the_geom)
    and a.id in (25758,24730,23875,25886,24442,26723);


texts:
Currently I'm using ArcMap>Image Analysis>Greyscale function to greyscale an OS VectorMap tiff. This turns hill or housing colours into greyshades as expected & shown in the image below.

Is there a way to greyscale the image but remove the grey fill on block shading like houses, hillsides or roads?


texts: Open WMTS query parameter encoding with "+" character in QGis

texts:
I am trying to create a WMTS endpoint with query parameters in the urls that contain a + character. So my resource url looks like this:

<ResourceURL format="image/png" template="http://localhost/wmts/{TileMatrix}/{TileCol}/{TileRow}.png?formula=(B8-B4)/(B8+B4)" resourceType="tile"/>

The formula query parameter contains a + character, which is not escaped when opening the endpoint in QGis. So in the html requests made by QGis, the formula is passed unaltered, and so the server interprets the + as white space. The specific calls look like this:

http://localhost/wmts/2/1/0.png?formula=(B8-B4)/(B8+B4)

But should look like this (plus sign escaped)

http://localhost/wmts/2/1/0.png?formula=(B8-B4)/(B8%2BB4)

To work around this, I tried to manually escape the character by changing the endpoint to

<ResourceURL format="image/png" template="http://localhost/wmts/{TileMatrix}/{TileCol}/{TileRow}.png?formula=(B8-B4)/(B8%2BB4)" resourceType="tile"/>

But in that case, QGis encodes the % symbol, leading to a dobule escaped formula that looks like this:

http://localhost/wmts/2/1/0.png?formula=B8%252BB4

So it looks like QGis escapes % but it does not escape + characters.

Is there any way to control this behavior?


texts: How we see all features of plugins in qgis 2.8.1?

texts: I have installed Qgis 2.8.1 but in plugins not installed, upgradable are not display. how we add new plugin?


texts: Writing code for Thematic map in Mapbasic?

texts:
I would like to write code for creating a thematic map.

The GUI command in Mapbasic window  shows many lines in order to set up the different colours which will be used on the different polygons (i.e. Brush (2, 16760, 2300)).

Is it possible to generate this part automatically provided the user gives some parameters ?

Part of the code would look like the following:

Declare Sub Main()
Declare Sub sThematic

Sub Main()
  Call sThematic
End Sub

Dim nStep As Integer
Dim sSymbole As Brush
Dim sSymbol(300) As Float
Dim counter As Integer
Dim sTableName As String
Dim sDataCol As Integer
Dim nBreak As String
Dim sStatement As String
Dim sStatementStart As String
Dim sStatementMid As String
Dim sStatementEnd As String

SUB sThematic

sSymbol(1) =   MakeBrush (49, 12632256, 12)
sSymbol(2) =   MakeBrush ("+sDefaultSymbol+", 16752895, 8)
sSymbol(3) =   MakeBrush ("+sDefaultSymbol+", 9445631, 8)
sSymbol(4) =   MakeBrush ("+sDefaultSymbol+", 255, 8)
sSymbol(5) =   MakeBrush ("+sDefaultSymbol+", 65280, 12)
sSymbol(6) =   MakeBrush ("+sDefaultSymbol+", 16776960, 12)
sSymbol(7) =   MakeBrush ("+sDefaultSymbol+", 16744448, 12)
sSymbol(8) =   MakeBrush ("+sDefaultSymbol+", 16711680, 12)

sStatementStart = "shade window " + winID + " " + sTableName + " with " + sDataCol + " ranges apply all use all Symbol (49,12632256,12)  "
sStatementMid = ""
sStatementEnd = " default Brush (2,8421504,16777215)  # use 2 round 1 inflect off Symbol (0,0,0) at 5 by 0 color 1 #"

nStep = 8
For counter = 1 to nStep
   sStatementMid = sStatementMid + nBreak(counter) + ": " + nBreak(counter + 1) + sSymbol(counter)
Next

sStatement = sStatementStart + sStatementMid + sStatementEnd

Run Command sStatement

End Sub


texts: Creating line from points and counting the number of points on that line

texts:
I have a shapefile that contains a number of points, each point represents a tree. I would like to calculate how many trees there are per row.
Hence my question (I think this is the best method), how do I automatically draw a line along a row (ideally of a baseline direction) and then calculate how many points are on that line (row)?


texts: many-to-many relate on two polygon layers in ArcMap?

texts:
I'm using ArcMap 10.3.  (Fairly new to GIS in general!)

I have two polygon shapefiles - school district and council districts - where the polygons overlap irregularly (meaning a school district can intersect with several council districts and vice versa). There is no key field to link them in any way.

I'd like to create a table that relates the two--really for analysis outside of ArcMap.  Someone asked this for QGIS, with some passing references to ArcMap here: Creating many-to-many table from intersection between two polygon layers in QGIS?

But I'm hoping for some direction on ArcMap.  Seems like it should be doable.


texts: Changing list-item css class using ng-class when mousing over Leaflet Map Markers?

texts:
I have an app with a map/markers on the right and a list-item menu on the left with info about the markers (like Yelp or Foursquare).

With some of my beginner hacking, I got the hover events to sort of work:

But it's odd, the list-item (pink background on hover) only works when I mouseout of the marker.  I'm trying to set it up so that when you mouse over the marker, the appropriate list-item's background changes and when you mouseout, it goes back to white.  Most of the other ng-class examples/questions I read through seem to work quite differently (they're going for a different functionality).

Ok, to the code:

HTML

The key part there is the ng-class="{hover: $index == hoveritem}"

Now I'll show you where hoveritem comes from

Controller

$scope.hoveritem = {};

    function pointMouseover(leafletEvent) {
        var layer = leafletEvent.target;
        //console.log(layer.feature.properties.id);
        layer.setIcon(mouseoverMarker);
        $scope.hoveritem = layer.feature.properties.id;
        console.log($scope.hoveritem);
    }

    function pointMouseout(leafletEvent) {
        var layer = leafletEvent.target;
        layer.setIcon(defaultMarker);
    }

    $scope.menuMouse = function(show){
        var layer = layers[show.properties.id];
        //console.log(layer);
        layer.setIcon(mouseoverMarker);
    }

    $scope.menuMouseout = function(show){
        var layer = layers[show.properties.id];
        layer.setIcon(defaultMarker);
    }

    // Get the countries geojson data from a JSON
    $http.get('/json/shows.geojson').success(function (data, status) {

        angular.extend($scope, {
            geojson: {
                data: data,
                onEachFeature: function (feature, layer) {
                    layer.bindPopup(feature.properties.artist);
                    layer.setIcon(defaultMarker);
                    layer.on({
                        mouseover: pointMouseover,
                        mouseout: pointMouseout
                    });
                    layers[feature.properties.id] = layer;
                    //console.log(layers);
                }
            }

        });
    });

}]);

So mousing over a marker
(layer.on({
                        mouseover: pointMouseover,
                        mouseout: pointMouseout
                    });)

fires the appropriate functions which then changes the icon colors.

I connected the layer.feature.properties.id; to $scope.hoveritem so that my HTML can then use that as the index (for ng-class).  When you mouseover a marker, it then feeds the marker id through to $scope.hoveritem which then it turn goes into the $index part of the HTML, thus changing it's CSS class.

But something is awry.  It only changes to the correct list item on mouseout instead of mouseover.  Furthermore, I can't figure out to get it to return to the default white state.  None of the list items should look active if the mouse is not on a marker.


texts: Matching land parcel polygon corners to corners on World Imagery basemap

texts:
I'm using ArcGIS Desktop 10.3.

I have a land parcel/cadastre layer of polygons which do not match up with the street corners on the World Imagery satellite basemap of ArcGIS.

Land Parcel layer is using, WGS_1984_Web_Mercator_Auxiliary_Sphere.

How do I get the Land Parcel corners to line up with the World Imagery basemap street corners?


texts:
I am just learning PostGIS and have a simple problem I can't get the answer for.

I want to select the geometries that intersect with given multiple geometries.

To make it simpler, let's say I want to select the line(s) which intersect with polygon1 and polygon2 in the image below.

I came up with something like this, but it's not working. What am I missing?

SELECT lines.id, lines.geom
FROM lines, polygons
WHERE(st_intersects(lines.geom,polygons.geom) and polygons.id = 1) and
     (st_intersects(lines.geom,polygons.geom) and polygons.id = 2)


texts:
I have two objectives to be maximised in ArcMap. One is wind speed and other is distance to residentials. I have the maps of wind speed and distance to residentials.
Is it possible to compare pixel 1 and pixel 2 and find if pixel 2 dominates or has better wind speed and distance to residentials?
If so pixel 2 will be eliminated. Likewise, I want to compare pixel 3 with pixel 4, pixel 5 with pixel 6 and all pixels in the study area.


texts: Looking for an Elko County Nevada Land Ownership Shapefile

texts:
I am looking for a Nevada Land Ownership shapefile.

County and state GIS don't seem to have anything beyond pdf.

Are there any public or commercial avenues to acquire land ownership in Nevada?

I'm really just needing public vs. private info.


texts:
I have a Tiff file that has no data value = -3.40282347e+38. I want to filter no data and get raw values with gdal. I tried with below code:

import numpy as np
from osgeo import gdal, gdal_array

dataset = gdal.Open("path/to/file.tif")
array = dataset.ReadAsArray()

lst = []
for v in array:
    if v < 0:
        lst.append(v)

print lst

but I received this error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/tmp/tmp3B0d5H.py", line 10, in <module>
    if v < 0:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

How can I do it?


texts: This is probably easy, but haven't found a solution. Convert Labels to Annotation is greyed out on right clicking a feature in the TOC. I have no reference scales set. Thanks.


texts: How to seperate class labels of hyperspectral images using R?

texts:
The image has different classes and different band, where the 220 bands are separated by the following code

   library(RStoolbox)
   library(raster)
   library(rgdal)
   library(caret)
   library(hyperSpec)

   img<-brick("19920612_AVIRIS_IndianPine_Site3.tif")

   names(img)
   names(img) <- c(paste0("B", 1:220, col = ""))
   plot(img$B1) ## plotting Band 1 image

How can we separate class labels such as alfalfa, corn-min etc from the above band?


texts: Problem or ERROR in converting from raster to vector!

texts:
I wonder if anyone could help me out, i am kind of a beginner at QGIS.
I´ve built a DEM for a project that i am developing than i´ve made the slope raster with the terrain analysis,it seemed that worked everything fine. But for further development of my work i ll need to vectorize some slope classes,but unfortunately every time i go Raster>Conversion>Polygonize; i get a error message (image with the error attached to this question).
I thought that could be something wrong with the DEM that i ve produced.Just to check it  i downloaded a DEM from ASTER DEM and  i tried to do same with these and the same error occurs!
Could be a software bug? Or am i doing something wrong?
Or there is another way of doing this?


texts: Spike and crater detection in DSM using machine learning / deep learning / computer vision techniques

texts:
I'd like to try machine learning / computer vision / deep learning approaches for detecting artifacts in Digital Surface Models, and am specifically interested in doing this for common artifacts like spikes and craters. So far, I haven't found very helpful research papers or datasets for this, and am looking for any pointers related to this.
Edit: I'm unable to directly share data or direct examples since the data is proprietary and I don't have permission to share it, and also because it's (largely) unlabeled. Some links that could provide more context/examples of what I'm trying to do:

Generating Spike-Free Digital Surface Models from LIDAR

"Errors in DEMs are usually classified as either sinks or peaks. A sink is an area surrounded by higher elevation values and is also referred to as a depression or pit. This is an area of internal drainage. Some of these may be natural, particularly in glacial or karst areas (Mark 1988), although many sinks are imperfections in the DEM. Likewise, a spike, or peak, is an area surrounded by cells of lower value. These are more commonly natural features and are less detrimental to the calculation of flow direction." - source


texts:
We are building a process which will allow users to submit a table of GPS data to update our central repository for all GPS data. Let's say that the target is simply a SQL data table. At first I thought about trying to accomplish this with an SQL INSERT wrapped in a Python script (exposed a GP tool). However, a colleague suggested that we somehow just get the table to the server and let integration services handle the Insert/append to the SQL table.

My question is then, how to create a simple file retrieval service? I'm comfortable with the Python/GP service aspect of ArcGIS server. Basically it would be just the user submitting a table and placing it in a known location on the server.

EDIT: I should add i'm not looking to upload a file from HTTP location; it must allow access to 'local' files on the user's system.


texts: Programatically take screenshot of map extent in ArcMap 9.3

texts: Is it possible to programatically capture a screenshot of the current map extent (the map area with any active displayed layers) in ArcMap with the ArcObjects .NET SDK?


texts: Downloading ArcGIS Online symbology to ArcGIS for Desktop?

texts:
I would like to transfer the symbology from my ArcGIS online maps to desktop. Each time I download it reverts to default.

Does anyone know of a way to do this, or perhaps there is a symbology package I need to download?


texts:
I've been trying to use the clipper tool in QGIS 2.18.0. However the following error message pops up;

Output driver `GMT' not recognised or does not support
direct output file creation. The following format drivers are configured
and support direct output:
GeoRaster: Oracle Spatial GeoRaster
VRT: Virtual Raster
GTiff: GeoTIFF
NITF: National Imagery Transmission Format
HFA: Erdas Imagine Images (.img)
ELAS: ELAS
MEM: In Memory Raster
BMP: MS Windows Device Independent Bitmap
PCIDSK: PCIDSK Database File
PCRaster: PCRaster Raster File
ILWIS: ILWIS Raster Map
SGI: SGI Image File Format 1.0
Leveller: Leveller heightfield
Terragen: Terragen heightfield
netCDF: Network Common Data Format
HDF4Image: HDF4 Dataset
ISIS2: USGS Astrogeology ISIS cube (Version 2)
ERS: ERMapper .ers Labelled
RMF: Raster Matrix Format
RST: Idrisi Raster A.1
INGR: Intergraph Raster
GSBG: Golden Software Binary Grid (.grd)
GS7BG: Golden Software 7 Binary Grid (.grd)
PDF: Geospatial PDF
MBTiles: MBTiles
MRF: Meta Raster Format
PNM: Portable Pixmap Format (netpbm)
PAux: PCI .aux Labelled
MFF: Vexcel MFF Raster
MFF2: Vexcel MFF2 (HKV) Raster
BT: VTP .bt (Binary Terrain) 1.3 Format
LAN: Erdas .LAN/.GIS
IDA: Image Data and Analysis
GTX: NOAA Vertical Datum .GTX
NTv2: NTv2 Datum Grid Shift
CTable2: CTable2 Datum Grid Shift
KRO: KOLOR Raw
ROI_PAC: ROI_PAC raster
ENVI: ENVI .hdr Labelled
EHdr: ESRI .hdr Labelled
ISCE: ISCE raster
ADRG: ARC Digitized Raster Graphics
SAGA: SAGA GIS Binary Grid (.sdat)
DB2ODBC: IBM DB2 Spatial Database
GPKG: GeoPackage

I'm saving all my data and layers in one folder.
My layers are all in the same coordinate system.
I also tried saving the points I want to interpolate as a shapefile first before clipping
I've been using mask as the clipping mode
I tried saving it as a .tif file and all the other recommended extensions in the error message but it still doesn't want to work.


texts: Re-initializing OpenLayers Map with different projection

texts:
Is it possible to re-initialize OpenLayers Map with a different projection ?

Scenario:

OL Map with projection 4326 is created, along with click, measurement, drawing handlers.
A button to switch to projection 3857 exists
When the button is clicked, we 'ideally' want to re-initialize the map to 3857,
without re-init all the handlers (click, measurement, drawing etc), because the tiles are are in 3857.

How do I do this?


texts: GeoServer raster tiles are individually contrast-balanced: how to turn this off?

texts:
I'm running GeoServer locally as a servlet (in JBoss), and I'm serving a GeoTIFF file through the WMS protocol.  In a webpage running an OpenLayers map, I have the GeoTIFF-WMS layer overlaid above an OpenStreetMaps base layer.  The problem is that it looks like this:

The GeoTIFF is clearly being tiled (good) and each tile is normalized to a given brightness level (bad).  When I zoom in, I can see a similar effect taking place on the smaller tiles, and I know that the original GeoTIFF doesn't have this blockiness: it's a diagonal strip of 16-bit greyscale surrounded by zeros that ought to be transparent.

How do I turn off this brightness-normalization feature (or what is it called, so that I can Google it)?

Even better, how do I tell GeoServer that pixels with zero values ought to be transparent?


texts: Set label precision with QgsRendererRangeLabelFormat class

texts:
I'm trying to control the number of decimal places represented for a layer in my project. I can do this in the symbology tab like this:

I want to do this in PyQGIS. Here is an example of how I would symbolize my layer:

uri = '/Users/ep9k/Desktop/SandraMonson/TestZips.shp'
join_layer = iface.addVectorLayer(uri, 'Patients by Zip Code', 'ogr')
target_field = 'Patient_Da'

def apply_graduated_symbology():

    myRangeList = []

    symbol = QgsSymbol.defaultSymbol(join_layer.geometryType())
    symbol.setColor(QColor("#f5c9c9"))
    myRange = QgsRendererRange(0, 3, symbol, '0-3')
    myRangeList.append(myRange)

    myRenderer = QgsGraduatedSymbolRenderer(target_field, myRangeList)
    myRenderer.setMode(QgsGraduatedSymbolRenderer.Custom)

    join_layer.setRenderer(myRenderer)

    print(f"Graduated color scheme applied")

apply_graduated_symbology()

The class QgsRendererRangeLabelFormat seems to be what I need to control the precision of these labels, as it has a method setPrecision.

What I do not understand is, it seems QgsRendererRange and QgsRendererRangeLabelFormat have no relationship to each other according to the documentation, the PyQGIS console, etc. In fact, QgsRenderRangeLabelFormat seems to have no parent or child classes and I can't figure out how this class relates to anything else.

It seems like I should be able to do something like this:

label = QgsRendererRangeLabelFormat('0-3', 2)   #2nd argument sets the precision
myRange = QgsRendererRange(0, 3, symbol, label)

But this really doesn't make sense because I create a string ('0-3') and then setPrecision, which should affect the precision of a float.

Here is the help from the QGIS console regarding arguments. Why is the first argument a Qstring?:

How do I use the class QgsRendererRangeLabelFormat?


texts: How to create equidistant points in QGIS along a polygon?

texts:
I am looking for a solution to distribute with a python script equidistant points in QGIS along a polygon. I tried to use the following code example:

http://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/

Using it in my example (running over each polygon from a polygon shape file):

# Reading a polygon layer
pp_layer = QgsVectorLayer("poly.shp" , "polylyr", "ogr")

# create layer to store result"
pntl = QgsVectorLayer("Point", "name", "memory")
pr = pntl.dataProvider()

for feature in pp_layer.selectedFeatures():
   geom = feature.geometry()

   currentdistance = self.distance
   feats = []
   # geom.asPolyline()
   while currentdistance < length:
      point = geom.interpolate(currentdistance)
      print point
      fet = QgsFeature()
      fet.setGeometry(point)
      feats.append(fet)
      currentdistance = currentdistance + self.distance
   pr.addFeatures(feats)

However, when running the code, if fails. I got an error message in the python console

Traceback (most recent call last):
File "<input>", line 1, in <module>
File "/home2/fun.py", line 92, in run_fun
fet.setGeometry(point)
TypeError: QgsFeature.setGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType'

an I got an errormessage from qgis

    ERROR 1: IllegalArgumentException: LinearIterator only supports lineal geometry components

Does anyone has an idea, if i have an error in my code or is it not possible to solve the question with this approach? I am using QGIS 2.0.1.


texts: How to find X,Y coordinates of a single point in ArcMap 10.1?

texts: I have a shapefile of a road layout (name: US-101). I drew a line across the road section which intersects the left and right edges of the road and pavement marking (which is also the part of shapefile US-101). Now, I want to find the X, Y coordinates (location in arcmap) of the intersection points of Us-101 and drawn line. The projected coordinate system of dataframe and shapefile is NAD_1983_StatePlane_California_III_FIPS_0403_Feet. This means that the coordinates should be in feet. If I use identify tool, it selects the whole shapefile and gives some location but I am not sure if that represents the location of the particular point I clicked on OR the location of shapefile. Please help me finding the coordinates. Refer to attached image.


texts: QGIS-3 Controlling PDF Export Image Quality and Compression

texts:
Is there a way of controlling what image compression or quality is used when exporting a PDF in layout composer?  ArcGIS users will be familiar with the options available under the "Format" tab in the ArcMap "Export Map" menu.

Is there a similar level of control available in QGIS 3.10?
I found a similar question but it seems dead:
Setting compression method when exporting composer from QGIS to PDF


texts: How to view a video using map-tip QGIS 2.16

texts:
In qgis 2.16 was improved map tip (changelog.qgis.org/en/qgis/version/2.16.0/#map-tip-improvements)
I did some tests to viewing a youtube video using:

<Iframe width = "420" height = "315"
src = "https://www.youtube.com/embed/XGSy3_Czz8k?autoplay=1">
</ Iframe>

(See screenshot)
in the map canvas appears only a black rectangle.

where am I wrong? how do you see the video?

Win 10 64 bit; OSGeo4W


texts: Making offset to polyline using ArcGIS for Desktop?

texts:
Is there a way to make an offset to a polyline in ArcGIS?
i.e. something like the command "Offset" from AutoCAD.

https://www.youtube.com/watch?v=9A8CyoAMnjU

I can't use the Buffer tool because I need to be able to calculate the length.

Also the "New polyline" has to retain the attributes of the original polyline.


texts:
I am trying to draw polyline on map, but line is not adding to the map,please help me determine where I am making a mistake.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>Simple Map</title>
    <link rel="stylesheet" href="https://js.arcgis.com/3.15/esri/css/esri.css">

     <button data-dojo-type="dijit/form/Button" id="drawPoly">Draw Polyline</button>


texts: How to do small expansion into LINESTRING ends?

texts:
The geometric problem is similar to this other one, that can be described as "small extrapolation of the line ends".
The correct solution must to expand/extrapolate d meters (ex. 2 meters) of line ends, of lines of any length.

Illustrating by a real life example
Here we need a  small expansion into line ends to glue lines of a road network.  The magic function not exists, but in some cases we can use ST_Scale to illustrate in the context.
I am testing this algorithm with no changes, except to include ST_Scale:
CREATE VIEW v AS
  SELECT city_id, (st_dump).path[1] as poly_id, (st_dump).geom
  FROM (
    SELECT city_id, ST_Dump(ST_Polygonize(geom))
    FROM (
       SELECT city_id,
              -- this works only for some straight lines
              ST_Union( ST_Scale(geom,1.01,1.01) ) as geom
       FROM road_network -- suppose all LINESTRINGs
       GROUP BY city_id
    ) mergedlines
    GROUP BY city_id
  ) polys
;

This wrong use of ST_Scale() funciton expanded in 1% the length (average ~2m) of the lines... Seems that works with straight lines, but, of course, curved lines will be distorted. And also I not need "proportional expansion" but only "constant line ends expansion" (the small 2 meters extrapolation of line ends).
I not see how to use the @SzieberthAdam's function of his solution, replacing ST_Scale  of the example. The @Jayden and @EoghanM solutions must be transformed into a generic function to be applied into LINESTRINGs.
PS: on this problem the ST_SnapToGrid is valid but only for very small correction. Not a good solution, because cause geometric distortions, and  the ST_SnapToGrid(geom,w) not solves all cases with w=2m (need bigger), and it not shows convergence with w (destroys/collapses the network).


texts: Using multiple colors for labelling that depend on values in QGIS?

texts:
I would like to do labelling like this (3 rows):

   Location_name
   Time
   People

My csv has 5 features

   Location_name      Time   AM/PM   People  A/K

Based on the value of "AM/PM", I want the "Time" has different colors; based on "A/K" values, I want the "People" had different colors.

Can someone help?


texts:
I am trying to calculate acos in ArcMap (shp) using Python. Here is the code:

arcpy.CalculateField_management(fc_Arcs, "Angle", math.acos(0.88), "PYTHON_9.3")

If I set the value i.e. 0.88 I am getting the output. However, I want to set the formula to calculate angle for all records in the shp.

I've tried to set the formula:

arcpy.CalculateField_management(fc_Arcs, "Angle", math.acos('!X!/!Radius!'), "PYTHON_9.3")

and I've got the error "TypeError: a float is required"

What I am doing wrong?


texts:
I'm working with a KML file generated by Google Location History takeout. Its structure is slightly different from a usual KML: each waypoint is recorded with 2 same-level tags (<when> and <gx:coord>):

<?xml version='1.0' encoding='UTF-8'?>
<kml xmlns='http://www.opengis.net/kml/2.2' xmlns:gx='http://www.google.com/kml/ext/2.2'>
 <Document>
  <Placemark>
   <open>1</open>
   <gx:Track>
    <altitudeMode>clampToGround</altitudeMode>
    <when>2011-03-20T17:33:11Z</when>
    <gx:coord>37.7605394 55.7806361 0</gx:coord>
    <when>2011-05-29T15:53:14Z</when>
    <gx:coord>4.5734683 45.0497192 0</gx:coord>
    <when>2011-06-09T16:45:50Z</when>
    <gx:coord>4.5736926 45.0501430 0</gx:coord>
   </gx:Track>
  </Placemark>
 </Document>
</kml>

On my understanding, readOGR() can't work with this kind of data. Are there any other solutions?

P.S. Requested code and output:

ogrListLayers("file.kml")
# [1] "Layer #0"
# attr(,"driver")
# [1] "KML"
# attr(,"nlayers")
# [1] 1

readOGR("file.kml","Layer #0")
# Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
# use_iconv = use_iconv,  :    Cannot open layer


texts: Graphically select a subset of Shapely objects in a Jupyter/iPython notebook

texts:
This question is being cross-posted from stackoverflow.

The interoperability between Shapely and Jupyter/iPython is nice. I can do cool stuff like create a bunch of geometric shapes and view them in the notebook:

some_nodes = [[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2]]
some_boxes = []
some_boxes.append([some_nodes[0], some_nodes[3], some_nodes[4], some_nodes[1]])
some_boxes.append([some_nodes[1], some_nodes[4], some_nodes[5], some_nodes[2]])

from shapely.geometry import MultiPolygon, Polygon
MultiPolygon([Polygon(box) for box in some_boxes])

...and Jupyter will show me this:

Now THAT is cool! It has been especially useful to me for quickly viewing and editing, for example, polygons that make up 2D finite element meshes.

Sadly, the images produced are just static SVG graphics; there's no interaction built in. It would be helpful to be able to select a subset of those objects in the image using this same graphical interface in iPython.

More specifically, I'd like to be able to create a list and add some of the displayed polygons to it by, for example, clicking/selecting them, or dragging a lasso/box around them, and also perhaps remove them when clicking a second time.

I've looked into trying to do this with matplotlib or javascript and although I've had some initial success it's probably the kind of project that is out of my depth at my current level of knowledge/skill.

Since Jupyter is a somewhat sprawling tool with lots of features that I am probably not aware of, I am wondering if there are already existing solutions out there for interaction of this sort in the context of a Jupyter notebook?

UPDATE: It looks like I'm going to have to create something myself. Happily, this tutorial is going to make that a LOT easier.


texts:
Just want to know if there is a tool to do a many-to-many join in QGIS?
I have X number of polygons that I want to join with a CSV table which contains Y number of entries. Essentially I need to have Y copies of feature X1, X2, ... ending up with X*Y features.
For now, I'm using a workaround which is to duplicate the entries in the CSV table X times and then do a one-to-many join.


texts:
I have imported a .tiff image of water pipelines(blue, red, purple), into ArcMap. I want to convert just the water lines into polyline vectors for measuring. Will Raster to Polyline work for this, or is there more I need to do?

Is there some sort of code to have ArcGIS recognize the colors and lines?


texts: QGIS not interacting with iface from python class

texts:
So I'm working on something where i need to use threads from python console of QGIS.
The following line of code adds layer successfully

iface.addVectorLayer("D:\ArcGIS Misc\Atlas\points.shp", 'Points', 'ogr')

But when i try to run the following code:

from threading import Thread
class AddLayer(Thread):
    def __init__(self, iface):
        Thread.__init__(self)
        self.iface = iface

    def run(self):
        self.iface.addVectorLayer("D:\ArcGIS Misc\Atlas\points.shp", 'Points', 'ogr')
        print 'added'

AddLayer(iface).start()

then i only get added printed on the console but the layer is never added. When i do it from a class which is not inherited from Thread then it works fine. Help Needed


texts: Make lines from points in Geometry Generator

texts:
I have a data set with a large number of start and end positions.  I am trying to generate lines between start and end positions using the Geometry Generator but it does not work.  I get no errors when running it, but the lines do not appear. I have found several answers to a similar qestion, e.g. this one, : Creating a line from two points in QGIS  I select geometry Type as "LineString.   Here is the code I have been trying to use.  make_line(make_point("shootlon" ,"shootlat"), make_point("haullon", "haullat"))

Small screenshot of the data, which is a text ("txt") file.
  What could be the problem?   Some of the older plugins e.g. "Points2one", do not work in QGIS 3.4.4.


texts: Convert a Tif raster file to vector shp file in QGIS

texts: So I've been trying to convert a georeferenced Tif raster file into a vector shp file using QGIS.  It's seems to complete this but it ends up giving me a feature class that has over 500,000 features, almost like it's creating a row per pixel.  So with that many features it's very hard to work with and when I do finally get it to Classify like into 8 classes of natural breaks, it doesn't look like the original tif file.  My end goal is to only have a few (less than 20) polygons based off the raster not thousands.  Could it be that the Tif file I'm using is already pre-processed or something?


texts: Is there a QGIS plugin that can transform tabular data to graph?

texts:
I have a CSV file with water level measurments in a well, taken once a year for 40 years.
Is there a way to show this data as a graph (x axis being the date and y axis the water level) in QGIS?


texts: Using ArcGIS sdk 10.2.2 for Android with gradle

texts:
Has anyone got the ArcGIS SDK from Esri to work with Android Studio(gradle)?. So far I have the done the following

Extracted the ArcGIS Android SDK(10.2.2).
Added the .so files (armeabi-v7a & x86) to the jniLibs folder as per this guide
http://ph0b.com/android-studio-gradle-and-ndk-integration/
Added arcgis-android-api.jar, arcgis-android-app-toolkit.jar, jackson-core-lgpl-1.9.5.jar, jackson-mapper-lgpl-1.9.5 and jcifs-1.3.17.jar to the libs folder

Here is my gradle file

buildscript {
repositories {
    mavenCentral()
}
dependencies {
    classpath 'com.android.tools.build:gradle:0.9.+'
}
}
apply plugin: 'android'

repositories {
mavenCentral()
}

android {
compileSdkVersion 19
buildToolsVersion "19.0.3"

defaultConfig {
    minSdkVersion 16
    targetSdkVersion 19
}

productFlavors {

    armv7 {
        ndk {
            abiFilter "armeabi-v7a"
        }
    }

    x86 {
        ndk {
            abiFilter "x86"
        }
    }

}

packagingOptions {

    exclude 'META-INF/LICENSE'
    exclude 'META-INF/LGPL2.1'
    exclude 'META-INF/NOTICE'
    exclude 'META-INF/license.txt'
    exclude 'META-INF/notice.txt'

}

}

A Screenshot of my Project Structure

The gradle file compiles successfully but when I  run the application it crashes with the following error

  Caused by: java.lang.NoSuchMethodError: no method with
  name='requestRender' signature='(Z)V' in class
  Lcom/esri/android/map/MapSurface

Am I missing something here?


texts: "...is not a function" error in Earth Engine - converting array to list

texts:
I am trying to convert the result of a region reduction by ee.Reducer.fixedHistogram() to an ee.List object so that I can .map() a function over it, but I get an error:

  [object].toList is not a function

The API doc for ee.Reducer.fixedHistogram() says that the output is an Nx2 ee.Array object. ee.Array objects have a .toList() method - why does it not work?

Here is my code:

// Define an image - select the SWIR1 band.
var img = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_038029_20180810')
  .select('B6');

// Define an area of interest.
var aoi = ee.Geometry.Polygon(
  [[[-110.5241853160436, 44.72231727808245],
    [-110.5241853160436, 44.44355785401221],
    [-110.06001295276235, 44.44355785401221],
    [-110.06001295276235, 44.72231727808245]]], null, false);

// Show image and AOI on the map.
Map.centerObject(img, 8);
Map.addLayer(img, {min: 100, max: 3500});
Map.addLayer(aoi);

// Calculate a histogram for the pixels in the API.
var swir1HistDict = img.reduceRegion({
  reducer: ee.Reducer.fixedHistogram(0, 5000, 250),
  geometry: aoi,
  scale: 30,
  maxPixels: 1e13
});

// A dictionary was returned, get the ee.Array results.
var swir1HistArray = swir1HistDict.get('B6');

// Convert the array to a list so that it can be mapped over.
var swir1HistList = swir1HistArray.toList();
print(swir1HistList);

Code Editor script


texts: Get the CRS of a raster file in a standalone script - Qgis Python

texts:
In a standalone Qgis python script I would like to get the CRS of a raster file. The idea is that I ultimately want to reproject it (probably using processing.runalg("gdalwarp ... ) but I don't want to hardcode the origin CRS in the code.

I have seen directions like this one but it all relies on the iface, which is apparently not straight-forward to use in a standaline script. Or this one, but it applies to shp and I didn't manage to adapt it to a raster layer

Is there an easier solution?


texts: Can I programmatically output the same map with different shp files?

texts: I have 39 different .shp files of the same locations but different data values and I need to output a map for each of them.  Is there a way to build a map template and then run a process that will loop through each .shp file and output a map, either .pdf or .png, of the map template with the individual .shp files?


texts: QGIS 3.2 Auxiliary Storage feature. Can it be bypassed or deactivated?

texts:
Up until QGIS 2.18, all our text placement and rotation values were set to point to the spatial layer table in MSSQL, where we added an ANG (angle) and X, Y field for rotation and placement value storage into our text layer tables.

After installing QGIS 3.2, it seems any label tool we use invokes a prompt to specify a primary key in order to set up fields in the new Auxiliary Storage feature.  This has undesirable results as all text is set to zero rotation, except the target label that was rotated.

We do not wish to use this new Auxiliary Storage feature and want to continue using the field values, from the base tables, to store our text placement and rotation data.

Attempting to revert the rotation source to the previous field in our base table causes QGIS 3.2 to crash.

Is there a way to turn off or bypass this new Auxiliary Storage feature?


texts:
I would like to display a date field as a label. The data is from a shapefile. The field in question is a date field.

When I use the date field for the label, it displays in the yyyy/mm/dd format. I would like to get it to the dd/mm/yyyy format.

Is this possible. I know in MapInfo I would use the FormatDate function, I wondered whether something similar existed in QGIS.


texts: Cloud masking for sentinel surface reflectance multi-images.But I cant apply for multi-images using Google Earth Engine JavaScript API

texts:
I want to do cloud masking for sentinel surface reflectance images but I cannot apply for multi-images,It was done for first image.
How could I do that using the Google Earth Engine JavaScript API?
var cld = require('users/fitoprincipe/geetools:cloud_masks')

var s2SR = ee.ImageCollection('COPERNICUS/S2_SR')
              //filter start and end date
             .filter(ee.Filter.calendarRange(2013,2021,'year'))
             .filter(ee.Filter.calendarRange(1,12,'month'))
             //filter according to drawn boundary
             .filterBounds(geometry)
             .filterMetadata('CLOUD_COVERAGE_ASSESSMENT', 'less_than',50);
//var first = ee.Image.constant(0) // init image

//var first = ee.Image(s2SR.get(0))
var test_image = s2SR.first()

Map.addLayer(test_image, {bands:['B4', 'B3', 'B2'], min:0, max:8000}, 'test image')

// adding mNDWI function
var addMNDWI = function(image) {
  var mndwi = ee.Image(image).normalizedDifference(['B3', 'B11']).rename('MNDWI');
  return ee.Image(image).addBands(mndwi);
};
 var mndwı=ee.ImageCollection(test_image)
             .map(addMNDWI);
 Map.addLayer(mndwı, {bands:['B4', 'B3', 'B2'], min:0, max:400}, 'mndwı')

var masked = cld.sclMask(['cloud_low', 'cloud_medium', 'cloud_high', 'shadow'])(test_image)
Map.addLayer(masked, {bands:['B4', 'B3', 'B2'], min:0, max:3000}, 'masked')


texts: Using QGIS Graphical Scale Bar Parameters

texts:
I am using QGIS 3.8.3.

Why in the composer layout, when I insert the graphic scale bar, both options (Cap style and Alignment) on the Item properties -> Display menu are disabled and the scale bar is always left aligned?


texts:
I have been making calls to a previously free API using esriRequest

var layersRequest = esriRequest({
    url: base_url,
    content: { "latitude": lat,
                "longitude": long,
                "riskCategory": risk,
                "siteClass": siteClass,
                "title":"Default"},
    dataType:'jsonp',
    handleAs: "json",
    callbackParamName: "callback"
}); ...

Now the API requires a key; which I have:

const key = 'xI234fg543example23';

The above is just an **example key; not real

The API documentation says to "Include your key in the api-key field of your request headers to authenticate like this:

curl --header "api-key: [your_api_key]" -X GET "https://api-hazards.atcouncil.org/public/v1/[load_category].json?lat=[lat]&lng=[lng]"

How do I correctly apply this using esriRequest?


texts: Calculate Distance between a Polygon and a Point Along a Path

texts:
I am very new to ArcGis and am trying to figure out how to do this.

Basically, this is what I am trying to do:

I have a map of Canada and in it I have polygons representing areas (each about the size of a municipality). I also have points describing the location of every airport in Canada and the network of roads all across Canada. What I need to do is find the distance from each area to the nearest airport using a road.

What I tried to do already was find the centroid of each polygon because I thought that I could use that to calculate the distance between two points along a road. But then I realized that the centroid is not necessarily connected to a road (especially in the northern parts of Canada) therefore I cannot use that.

I would really appreciate some sort of guidance as to how I would accomplish this as I am very new to ArcGIS and have not found much to help me (without confusing me) online.


texts: Euclidean distance process is instantaneous without mask but runs for hours when boundary mask is added

texts:
I have shapefiles of lines and points which I need to create a Euclidean distance raster of using the Euclidean distance tool in ArcMap.
The size of the area is about 37,000 km². The boundary shapefile (polygon) is 35 mb in size. I set the cell size to 500 m and the maximum distance of calculation to 80,000 m and the process run time is a few seconds.
However, when I add the region boundary as a mask, the process takes over 2 hours to complete. I have run this task before for another area of comparable size, using the boundary mask of that area and the process ran for only a few seconds.
Why is this taking so long with my area now?
Could there be an issue with the boundary shapefile?
I have already repaired geometries.


texts:
I'm trying to upload a CSV-file with coordinates of about 80,000 lines, when ever I uploading the file to QGIS, for some reason, the coordinates don't display right and spreading equally. But when I upload one by one (for testing) the position is correct.

It does not matter from me the map I just want use intersection with polygons and points in order to get values for points.


texts: How to use raster calculator on WMS layers?

texts:
I'm using lidar data (UK Environment Agency) from a WMS server and would like to use this in the raster calculator. When I go to the raster calculator the WMS layers aren't an option.

The data is also available for download in ASC format, which can be made into virtual raster format that seems to work with the raster calculator. However, the lidar dataset can only be downloaded this way in small sections and it would be very time consuming to do this.

Is there a way to either:

1 - get the WMS layer to work in the raster calculator?

2 - download all the data from the WMS server and save it in a format that works with the raster calculator?


texts: How to add marker while click on geosearch location?

texts:
I am using the archgis map and basemap street.1) Fist i click on the geosearch button and then show the search location when i click on the search location one function fire which is shown below:

function chooseAddr(latt,longt)

{
$('#linkid').attr('data-latt', '');
var combine=latt+','+longt;
$('#linkid').attr('data-latt', combine);
 $('#linkid').click();
 //abc();
}

This is written on index.html page and on js file

$("#linkid").on("click",    function() {
    alert('inside the abc');

    var dataLatt = $('#linkid')=='undefined'||$('#linkid')==''?'0,0':$('#linkid').data('latt');
    var splitdata = dataLatt.split(",");
    var latt = splitdata[0];
    var longt = splitdata[1];
    addmarker(latt, longt);

    //addmarker('28.833333', '76.666667');

});
function addmarker(latt, longtt) {

    polygraphlayer.clear();

    var pt1 = new Point(Number(longtt), Number(latt));

    var markersymbol1 = new PictureMarkerSymbol('images/marker.png', 10, 15);

    polygraphlayer.add( Graphic(pt1, markersymbol1));

     //map.add(polygraphlayer);

    //map.graphics.add(polygraphlayer);

}

Its show the marker first time next time when i click on the search with different lat, long marker is not shown with same functions


texts: Comparing points (latitude/longitude) w/address data (zip) to zipcode boundaries shapefile using QGIS?

texts:
I have a bunch of 'points' in a csv that contain a zip code and lat/long coordinates. I also downloaded the ZCTA boundary data from the census website. I successfully imported the point layer on top of the zip boundary layer, however, I'd like to determine if there are any 'points' that lie outside of it's associated zip boundary layer.

1) Points = Zip Code, Latitude, Longitude
2) Zip Boundary Shapefiles

Basically I want to determine which points' latitude/longitude (and zip code) does not line up with its corresponding zip boundary layer. I definitely know this is possible, and I can also do it manually by hand, but

How do I do this for a set of points using QGIS?


texts: QGIS layer in .svg or with a better resolution?

texts:
I have a georeferenced layer, built from a nice map drawn with Inkscape. But when I convert it to *.png, *.jpg or *.tif, there is a pixelisation of the map. After georeferencing, it gives a pixelised map. I tried different definitions, sizes, etc but it always has much pixelisation.

How could I import the same map, but with nice contours, like in a *.svg file? I have thought about converting it into a shapefile, but I haven't found out if it is possible?


texts: Cannot change zoom level parameter dynamically within gdal2tiles call in Python?

texts:
I am trying to do the following but it is falling over i.e. (has no error message
and doesnt start the tiling process). I think the parsing of the zoom level argument is falling over. I havent been able to parse the zoom level

NB. the vrtName variable contains the zoomLevel as string "0-5", "5-10", or "10-18" which was taken from the vrt names 0-5.vrt, 5-10.vrt and 10-18.vrt

data = vrtName.split("-")

minZoomLevel = data[0]

maxZoomLevel = data[1]

print("gdal2tiles tiling process started.....")

os.system('python "C:\Program Files\GDAL\gdal2tiles.py" -p geodetic -z ' + str(minZoomLevel) + '-' + str(maxZoomLevel) + ' C:/Marine/mosaic_tif/' + vrtName + '.tif ' + outputFolder + clientName + str(n)) n = n + 1

Is it possible to change the zoom level parameter dynamically in a gdal2tiles call within python?


texts: Efficient way to use a CSV to Define a Feature Class

texts:
I was provided an Excel Spreadsheet which had the definitions of fields for a feature class.  There were several worksheets in the file, each containing field definitions for different feature classes.  After a bit of cleanup, I converted each worksheet to a CSV with the following columns:

Field_Name
Field_Type
Field_Length
Field_Alias

Each of the CSV files had between 300 and 400 rows, which means there are that many fields (I know it's a lot of fields, but it's not my schema design).

I put together a bit of code to process the CSV files and create fields for the feature classes.  Note, the feature classes already exist in a geodatabase and the CSV files have the same name as the feature classes.

The following code works, but it takes over 20 minutes to create 3 feature classes with 300-400 fields each.

count = 0
for csv_file in CSVList:
    print csv_file
    print count
    with open(csv_file, 'rb') as csvfile:
        reader = csv.reader(csvfile, delimiter = ',')
        for line in reader:
            field_name = line[0]
            field_type = line[1]
            field_length = line[2]
            field_alias = line[3]
            arcpy.AddField_management(FCList[count],field_name, field_type, "","",field_length,field_alias)
    count+=1

Is there a more efficient way to create a feature class from a CSV file with field definitions?  Or, is arcpy.AddField_management as good as it gets for this sort of thing?


texts: Correct format of query in geoprocessing tool

texts:
In the code below, I keep getting the error An invalid SQL statement was used. [SELECT OBJECTID FROM LARA_CDC WHERE County IN  ['44', '50', '74', '76']].

I've tried several different ways of formatting the query, with no luck. The end goal is to create a new file geodatabase table from selected records of an existing one. Any ideas where I'm going wrong?

new_table = "Childcare"
new_table_path = os.path.join(table_locations, new_table)
fld = "County"
query = "{0} IN ".format(arcpy.AddFieldDelimiters(table_locations, fld))
query += " ['44', '50', '74', '76']"
....
arcpy.TableToTable_conversion(table_path, table_locations, new_table, query)


texts: Debugging Error: Algorithm gdal:dissolve not found in PyQGIS?

texts:
I have difficulties with the processing package. I try to run a simple QGIS dissolve in pyQGIS 3. When I run the line :

from qgis.core import QgsVectorLayer

import processing
from processing.core.Processing import Processing
Processing.initialize()

layer = QgsVectorLayer(layerPath, 'layer 1', 'ogr')

parameters = {'INPUT': layer,
              'OUTPUT': 'memory:'}
dissolved = processing.run('gdal:dissolve', parameters)

In QGIS Python Console, it works fine but when I try it in my script, it sends me the error message in the title.
I run my scripts with a batch that imports all necessary libraries.

The full message is following :

File "C:\Program Files\QGIS 3.2\apps\qgis\python\plugins\processing\tools\general.py", line 96, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:\Program Files\QGIS 3.2\apps\qgis\python\plugins\processing\core\Processing.py", line 125, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Error: Algorithm gdal:dissolve not found


texts: How can I plot the availability of (cloud-free) data in a satellite image collection in Google Earth Engine?

texts:
I have often worked with time series, including the masking of clouds in GEE, including using the ui.Chart.image.series() function to plot time series for the values of a given pixel. What I haven't managed to achieve yet, though, is not the plotting of the time series of some value, but of the actual time stamps available.
Here's an example framework:
// Function to remove cloud and snow pixels
function maskCloudAndShadows(image) {
   var scl = image.select('SCL');
   var mask = ((scl.neq(3)).and(scl.lt(7)));
   return image.updateMask(mask);
}

var startDate = '2019-01-01'
var endDate = '2019-12-31'

// Use Sentinel-2 L2A data - which has better cloud masking
var s2_collection = ee.ImageCollection('COPERNICUS/S2_SR')
    .filterDate(startDate, endDate)
    .map(maskCloudAndShadows)
    .filter(ee.Filter.bounds(testPoint)) //testPoint can be any point

var s2_chart = ui.Chart.image.series({
    imageCollection: s2_collection.select('B8'),
    region: testPoint
    }).setOptions({
      interpolateNulls: true,
      lineWidth: 1,
      pointSize: 3,
      title: 'Band 8 over time at test point',
      vAxis: {title: 'Band 8'},
      hAxis: {title: 'Date', format: 'YYYY-MMM', gridlines: {count: 12}}
    })
print(s2_chart)

If you select some point on the map, this script will produce a nice time series of the reflectance observed in Band 8. But instead of this, I would like to get a plot that indicates those time stamps, from which the plot is created, e.g. by a point or bar plot of the same value for every date that comes with a non-masked value of Band 8 (and maybe also a point or bar plot of different color for every masked value of Band 8). It would even be OK to export a CSV file containing just the dates of the all values that would have went into the time series plot. How can this be achieved?


texts: How do I run gdalwarp on a number of csv of coordinates?

texts: I have a CSV of Ground Control Points in NAD83 UTM 12N and I need to convert them to WSG84. Is it possible to do this with gdalwarp or something similar?


texts: Standardized format for contributing street closure data to navigation apps

texts:
I am wanting to put together an application that will digest georeferenced road closure data. I am wondering if any such format currently exists for notifying google maps or other navigation applications so that I am not recreating the wheel.

If there is not a format that is currently used would some form of GeoJSON be the best option?