texts: Reading data from raster layer obtained via WMS in QGIS

texts:
I need to use some geological data from the WMS service http://mapas.igme.es/gis/services/oneGeology/IGME_EN_Geology/MapServer/WMSServer using QGIS. I've connected to the service and I can see the provided layers and if I click in the map for information I obtain the data OBJECTID, SHAPE, etc. But now I need to do some processing using Python and I'm a bit confused.

I've selected the active layer via

lyr = self.iface.activeLayer()

and I've discovered that the layer is raster. This is confusing for me because when I ask for information of an element by clickink in the map I obtain the field SHAPE as Polygon.

Well, I accept the layer is raster, so I will try to read the data. If I ask for information as

print(lyr.extent())

I obtain the rectangle limits that enclose the data in the layer. But if I try

print(lyr.width(),lyr.height())

I obtain

(0, 0)

Is this output correct or should be the number of columns and rows of the raster? Also, if I try

ident=lyr.dataProvider().identify(QgsPoint(-8.3,43.0),QgsRaster.IdentifyFormatValue)
print ident.results()

I obtain

{}

but if I click in the map at coordinates (-8.3,43.0) I obtain correct results.

Is it possible to iterate over all pixels in a raster layer from a WMS service using Python in QGIS?


texts: I would like to use qgis similarly as a chart plotter for dynamic data that are created by a Fortran code (being basically a ship simulator.) I have two questions: how do I integrate Fortran into qgis the best way and how can I dynamically update the plot of the position?


texts: How to write feature types with different geometries using Schema(From Table) to FME FileGDB using FME?

texts:
I need to generate a geodatabase with 25 feature classes according to varying schema (geometry included) that are defined in an excel spreadsheet. My approach so far has to been to use the Schema Reader(From Table) and the FileGDB dynamic writer. However, when I run the workspace it returns a single feature class, despite the fact that I have dynamic schema enabled and specified to define feature classes based on the input table.

EDIT: I used the parameters as shown below and was returned unique geodatabase tables, however they did not carry the geometry specified in the input schema.

Is there an additional transformer needed to run this type of translation? Obviously I do not want to fan out by feature type because that only allows one geometry type.

See below for sample of input schema:


texts: How to center a map in Print Composer to the same center as in map window?

texts:
I can't figure out how can I center the map in Print Composer to layer/object/map center in the same way as in map window in QGIS 1.8. Let's say I have layer which I want have in the center of the map. In map window it is pretty easy - I just click on
"Zoom to layer extent".

Unfortunately when I create map frame in Map composer, the layer is
centered on width (good) but not centered on height (I guess it is aligned
to bottom).

If I use button Set to map canvas extent in Print Composer / Item
Properties / Extents it is again centered on width but this time aligned to
top.

I think in both cases it would be much better to have it centered. Is there
any way to have map centered in Print Composer to the same center as map in
map window?


texts: Update field's rows with a value of dictionary

texts:
I'm new in programming and I'm struggling with my code. I have a dictionnary of districts and communities for each, like this:

MyDict = defaultdict(<type 'set'>, {u'DISTRICT': set([u'community1', u'community2']),
u'DISTRICT2': set([u'community3', u'community6'])})

Now I want to add the list of values in a shapefile of districts, in the row corresponding to the right district: so the field in this shapefile called 'District' has to be the same as the key of my dictionary if I understood well.
Expected result in shp:

<District>   |   <Communities>              |
DISTRICT     |   community1, community2     |
DISTRCIT2    |   community3, community6     |

My code isn't working, I get no error message but nothing written in the expected field:

Communities = ['Communities']           #field to be written with values
District = ['District']                 # field with district in the shp

# Update in the communities field of the Districts shapefile
with arcpy.da.UpdateCursor(MyDistrictShapefile, Communities) as cursor1:
for newRow in cursor1:

  # Define a loop pour all the district name in the shapefile
  districtsLoop = [row[0] for row in arcpy.da.SearchCursor(MyDistrictShapefile, District)]

  # Looking for the same district as the key and update the communities field
  if str(districtsLoop[0]) == str(myDict.keys()) :
     newRow[0] = myDict[str(districtsLoop[0])]
  cursor1.updateRow(newRow)

Anyone has an idea of how I can solve my problem? It's my first code and I don't see how to fix this.


texts: I am trying to create a DEM from topographic contours vector. I have already sucessfully created a raster using the interpolation plugin. Firstly, I don't know why or how to ccontrol the output settings when creating a raster so that it can be  FLOAT32 compliant. Can anyone shed some light on this issue?


texts: Can anyone tell how Sumo's netconvert is using netoffset (lat,lon offset) when we convert osm file to .net.xml file in order to import in sumo ? Netoffset is present in .net.xml file in location tag.


texts:
I am trying to create a UTFgrid JSON with a Raster in PostGIS. The goal is to provide weather data for web maps.

I found some libraries:

https://github.com/Ecotrust/create-utfgrids

https://github.com/springmeyer/utfgrid-example-writers

Is there one way to create the UTGrid from PostGIS using Javascript?


texts: Reprojecting and retiling a large set of tiles

texts:
I am trying to understand the process to re-project and re-tile a set of raster (GeoTiff) but can't wrap my head around how to use the different GDAL tools to achieve my goal.

More concretely, I would like to transform the full set of SRTM tiles that are in the geographic projection (5x5deg) into a QSC (Quadrilateralized Spherical Cube) projection to have a set of 256x256px tiles for each cube face.

Do I need to first use gdal_merge to create one big GeoTiff that I can then reproject and tile using gdal_wrap? Is it possible to perform this directly from the source tiles without merging (seems wasteful)?


texts: How to create a 'dummy' DATETIME field for QGIS timemanager?

texts:
We have several 'time-aware' datasets, but the dates are stored as YEAR (2017), SCHOOLYEAR (20172018), etc. rather than a true DATETIME value.

In order to animate these datasets, is it possible to create a dummy DATETIME value using a virtual field in QGIS?

I've tried so far to create a virtual field using the following expression:

to_datetime(LEFT("SCHOOLYEAR",4) || '-' || '10' || '-' || '01')

Which looks correct:

...though in the attribute table there is an added T - which might be acceptable based on the formats I get in the error below...

But when I apply the layer to the timemanager, I get this error:

Any thoughts on what format / expression I should be using given my inputs?


texts: Splitting polygons with line grid in QGIS?

texts:
I have a polygon vector layer and a line vector layer in a form of grid (they are both in the .shp format). I need to split polygons by these lines.

I tried digitizing tools but the problem is that it deals only with 1 line at a time and I have more than 52,000 lines. It will take a ridiculously long time to split it manually. How can it be done for all the lines?

I have QGIS 2.18.


texts:
I'm a new user of ArcMap. And I'm not good at programming language.

There are the following shapefiles in workspace.

L01-1983_01.shp
L01-1983_02.shp
...
L01-1983_46.shp
L01-1983_47.shp
L01-1984_01.shp
L01-1984_02.shp
...
L01-1984_46.shp
L01-1984_47.shp
...
L01-2012_01.shp
L01-2012_02.shp
...
L01-2012_46.shp
L01-2012_47.shp
L01-2013_01.shp
L01-2013_02.shp
...
L01-2013_46.shp
L01-2013_47.shp

File name has the following rule.

L01-(year)_(city code).shp

I want to merge all shapefiles for each year, for example,

L01-1983_01.shp
L01-1983_02.shp
...
L01-1983_46.shp
L01-1983_47.shp
-----merge----    L01_1983.shp
...
L01-2013_01.shp
L01-2013_02.shp
...
L01-2013_46.shp
L01-2013_47.shp
-----merge----    L01_2013.shp

Below is my script.

list_year=["list_%d" % x for x in range(1983,2014)]
year=1983
for i in list_year:
    i=ListFeatureClasses("*year*")
    year=year+1
years=range(1983,2014)
for i in years:
    Merge_management(list_i,"L01_i.shp")

But I can't get merged files.

How can I write python script to merge shapefiles using a loop?


texts:
I would like to know what is the sense of using st_buffer operation. In PostGIS reference guide there is only the syntax and a short description of this operation.

Could you give some examples to understand the advantages of this operation?


texts: Detect and measure loop at end of linestring

texts:
In PostGIS, is there a good way of finding out whether a linestring loops back on itself and, in that case, how long the loop is?

I am rendering 1:50000 topographic maps from OpenStreetMap data, trying to emulate the style of the maps produced by the Swedish geodetic survey, Lantmäteriverket. On these, a turning circle, as often found at the end of forestry roads, is represented by a small white disc. In OSM data, however, editors simply often draw a loop at the end of the road (because that is often what it looks like). So, I want to detect a small loop at the end of a linestring, so that I can render it properly.

Since I am working with OSM data, I have topological information as well as geometric, so the "whether" bit is easy. Just see if the end node of the way is also an interior node. But how do I "chop off" the loop to determine its length (or diameter, or any other dimension)? There is ST_Line_Substring, but that requires proportional distance along the way, which I do not have. I do have node indices, so I could do ST_LineLocatePoint(ST_PointN(...)) but since I will always have two coincident points (start and end of loop) I doubt this will work.


texts: Python: cfitsio library not loaded when importing OSGeo library or using GDAL on Mac

texts:
I am trying to use the OSGeo library and GDAL Python bindings in a script for Python 3.7.7 on MacOS 10.15.4 Catalina. However, when I try to compile the script, I receive the following error as soon as I import from the OSGeo module:

Traceback (most recent call last):
  File "/usr/local/lib/python3.7/site-packages/osgeo/__init__.py", line 18, in swig_import_helper
    return importlib.import_module(mname)
  File "/usr/local/Cellar/python/3.7.7/Frameworks/Python.framework/Versions/3.7/lib/python3.7/importlib/__init__.py", line 127, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
  File "<frozen importlib._bootstrap>", line 1006, in _gcd_import
  File "<frozen importlib._bootstrap>", line 983, in _find_and_load
  File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 670, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 583, in module_from_spec
  File "<frozen importlib._bootstrap_external>", line 1043, in create_module
  File "<frozen importlib._bootstrap>", line 219, in _call_with_frames_removed
ImportError: dlopen(/usr/local/lib/python3.7/site-packages/osgeo/_gdal.cpython-37m-darwin.so, 2): Library not loaded: /usr/local/opt/cfitsio/lib/libcfitsio.8.dylib
  Referenced from: /usr/local/opt/osgeo-gdal/lib/libgdal.27.dylib
  Reason: image not found

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "test.py", line 1, in <module>
    import countries
  File "/Users/adamsamec/Dropbox/Work/Eurofly/eurofly-reverse-geocoder/countries.py", line 4, in <module>
    from osgeo import ogr
  File "/usr/local/lib/python3.7/site-packages/osgeo/__init__.py", line 32, in <module>
    _gdal = swig_import_helper()
  File "/usr/local/lib/python3.7/site-packages/osgeo/__init__.py", line 31, in swig_import_helper
    return importlib.import_module('_gdal')
  File "/usr/local/Cellar/python/3.7.7/Frameworks/Python.framework/Versions/3.7/lib/python3.7/importlib/__init__.py", line 127, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
ModuleNotFoundError: No module named '_gdal'

Here are the steps I've taken to install the osgeo and GDAL libraries, partially using the Homebrew package manager for Mac:

Downloaded and installed this GDAL 2.4 Complete package.
In Terminal, run the following commands:

pip3 install gdal
brew tap osgeo/osgeo4mac
brew install osgeo-gdal-python

Looks like the GDAL library is not properly installed because when I run gdalinfo, I get the following error:

dyld: Library not loaded: /usr/local/opt/cfitsio/lib/libcfitsio.8.dylib
  Referenced from: /usr/local/bin/gdalinfo
  Reason: image not found
zsh: abort      gdalinfo

Do you have any ideas what I am doing wrong?


texts: Automatic contrast stretch for each atlas element?

texts:
I'm working with QGIS 2.18 (windows) and the useful atlas tools. My covering
layer is a polygon layer (20 features). My atlas template is composed by 4
map elements representing different images raster (3 multibands) located on the same area.

I use the  option proposed in the tab called "object property" to lock the layer and the style of the layer (lock layers / lock style for layers). For each
raster, I want to be able to automatically apply and set the raster style
and enhance the contrast "local histogram stretch".

I join a screen shot representing my atlas template.

Is there a way to automatically assign a contrast stretch ("local histogram stretch") for all my map elements and for each of my 20 atlas sheet?


texts: Using raster file and Excel data for analysis QGIS

texts:
I'm using the population density raster from SEDAC (person/square km), I've clipped it to my area of interest (Vietnam) as shown in the picture below:

In addition to this, I have a vector layer of the provinces of Vietnam:

For each province I have a unique id and the name of the province.
Lastly, I have an Excel spreadsheet containing for each province, the waste production rate per person. I'd like to group the population data for each province, a bit like "join attributes by location" with a vector, but I'm unsure I can do this with raster files. I'm starting to use raster info; I've looked around but haven't found anything helpful.


texts: How to find nearest city using Geonames features codes

texts:
I've got about 300,000 lat/lng coordinates of events that I'm combining into a large map. The coordinates come from all over the world.

For each point, I want to find the nearest city of at least XYZ residents. I've written a Node script that uses the geoNames database of world cities and does all the calculations.

It all works fine but for the fact that there seems to be a far amount of overlap in some of these cities. For example, cities15000.txt lists both "North Atlanta" and "Atlanta," and there is both a "Manhattan" and a "New York City."

I was hoping I could limit the search to the highest level in the city hierarchy using geoNames' feature codes, but the documentation is not terribly clear. North Atlanta is "PPL" (populated place), Atlanta is PPLA (seat of a first-order administrative division). I thought I could limit to PPLA, but Baltimore is PPLA2 (seat of a second-order administrative division). There does not appear to be anything near Baltimore (say, "City of Baltimore" or something) that I'm missing.

Where do these codes come from? Is there a way to eliminate subsections of cities? I also don't mind using a different source, though I do find the geoNames data clean and easy to manipulate.


texts:
how can i add a layer upload panel tool in my suite sdk?
The main problem is when i click "upload layers" button. Appear errors in the firebug. I don't
know the reason why those errors appear

TypeError:b[e] is not a constructor...
TypeError:this.dataStore is undefined:this.dataStore.emptyText = this.dataStoreEmptyText; (in LayerUploadPanel.js line 327).

i'm not sure if i haven't added the necessary code
to app.js.
i want to upload mi own shapes from layeruploadpanel, but i don't how to do it. Anyone could help
me please?


texts: Determine use of pipes in a geometric network

texts:
I am trying to understand how to use the Geometric Networks in ArcGis 10.1

For the purpose of my experiment, I built a simple network composed of a few connected pipes, one source (left triangle) and one sink (right square)

Following Esri's help documents, I managed to create the geometric network dataset and set my source and sink features as such.

My final objective is to find the segment that is the most likely to experience overcharge. Intuitively, it appears that it should be the one linking the two "blocks" with the directional arrow.

Is there a way to measure the charge over each edge of the network if, for example, I want 100 liters of water to flow from my source toward my sink? (given that all the pipes are identical in capacity)

By the way, I am also wondering why the "Set flow direction" command was able to give the direction of only one pipe, given the position of the source and the sink.

Thanks fo your help!

EDIT for clarity:

If I take another simplified network as example, I would like to quantify the charge/pressure on eache edge. In this network, I inject 100 units in the network from the left, and it moves to the right. While the edges in the loops only support 50, the segment in the middle has to resist to 100 units. Can ArcGis model the load of edges in a network?


texts: Loading .sid files in QGIS 2.18.3 for Ubuntu

texts:
I just installed QGIS 2.18.3 on Ubuntu 14.04 and just realized that .SID files are not accepted as a valid file format for raster images.

Is it possible that there is a library missing from the distro from ubuntugis?


texts: Any restrictions on calling ArcPy from non-Esri web service?

texts:
I recently saw some success calling a python geoprocessing script from Node.js.  I was also thinking of using flask instead of an arcgis server geoprocessing service.

Does anybody know if there are restrictions on calling arcpy from non-esri web service from a license perspective?

Disclaimer: No advice will be interpreted as legal counsel, just looking for general knowledge.


texts: Can you use animated images as symbols in Mapinfo?

texts:
How is it possible to use image that blink or change (animate) as
symbol in mapinfo? I tried to create custom symbol that blink in
mapinfo but after I added it as symbol it didn't blink. Have I done
something wrong?


texts:
I´m working on a GIS-Project in the Defense & Intelligence sector.
The goal was to develop/create a GIS-Service that:

stores and manage GIS-Data (Map data, Height data.. )
has a repository for GIS-Data ( for versioning, etc.)
Synchronization of GIS-Data between Workstation/PC and Tablet (over Network + local with usb sticks/ hdd’s)
GIS Functionalities (Analysis, etc.) as a part of the GIS-service and as a local library accessible from different apps on the Tablet
OpenSource
C++ / Qt – API
Placing symbols on the map

So the main goal is to store, manage and share the maps between different apps on tablets, servers, workstations and pc’s.

Analysis:

Line of Sight
Routing
Terrain analysis

I already compared the most popular OpenSource Gis-Tools but the main issue here is that I have to use C++ and Qt.

I also put a lot of effort into trying to build an standalone application based on QGIS which led to a lot of problems like compiler Errors and other building failures.

So is there a EASY to use OpenSource GIS-Tool for developing a GIS-Application or do I have to spend money on professional GIS-Tools ?


texts:
I have multiple polygons w/ attributes associated with each of them (as shown in the image).  I need to bring these into Google Earth and have the attributes look something like they are shown in the other image, on them in google earth.  I have done this in the past and had no problem but currently when I try to save my polygons as KML and bring them into google earth I cannot click on them to bring up the attributes.

I am using QGIS.


texts:
I'd like to compare LAI (Leaf Area Index) data from the same region over different times using ui.Chart.image.series in Google Earth Engine.
I'm trying to set the axis to the same range, in order to compare the charts by using min Value and maxValue like this:
var chart2 =
    ui.Chart.image
        .series({
          imageCollection: lai_august_19_broadleaf,
          region: aoi,
          reducer: ee.Reducer.mean(),
          scale: 500,
          xProperty: "system:index"
        })
        chart2.setOptions({
          title: 'Mean LAI France broadleaf forest forest June-August 2019',
          hAxis: {
            title: 'Time',
            minValue: 0,
            maxValue: 6
          },
          vAxis: {
            title: 'LAI',
          },
        });

The chart however does not adjust the axis using that method. Is this the correct way or is there another one to adjust the axis?


texts: Using Select statement inside of array constructor in PostgresSQL?

texts:
I'm working with point coordinates in PostGIS and I have a problem. My goal is to generate linestrings by using ST_MakeLine geometry constructor (see snippet 1). Source data is pre-generated ST_MakePoint statements and can contain multiple points (see example from snippet 2). Due that, I have to use ARRAY constructor inside ST_MakeLine geometry constructor.

-- snippet 1 (example of geometry contstructor I'm trying to use)

SELECT ST_SetSRID(ST_MakeLine(ARRAY[<point_data_here>]),4326)

-- snippet 2 (data in column, type text)

ST_MakePoint(62.577460, 29.789088), ST_MakePoint(62.577431, 29.788905), ST_MakePoint(62.577109, 29.789399), ST_MakePoint(62.576793, 29.790482)

-- snippet 3 (this works well, returns geometry as expected)

SELECT ST_SetSRID(ST_MakeLine(ARRAY[ST_MakePoint(62.577460, 29.789088), ST_MakePoint(62.577431, 29.788905), ST_MakePoint(62.577109, 29.789399), ST_MakePoint(62.576793, 29.790482)]),4326)

My exact question is, how can I embed text data shown in snippet 2 in ARRAY constructor shown in snippet 1? If I do it by manyally copying and pasteing (see snippet 3), there are no problems. But how can I do SELECT statement (data from column which have type of text) inside that ARRAY constructor?


texts: Trying my first map by uploading a SHP file.  The instructions tell me that I can upload a SHP file, but I get an error that they are not allowed.  This happens whether I upload by browsing a local directory, or by dragging a file into the CartoDB window.


texts: Automatically creating 10,000 sets of random points for further analysis in QGIS?

texts:
I'm using QGIS and would like to know if I can automate generating multiple sets of random points and mean distances?

I have a data set consisting over 600 coordinate points for Mistletoe across Greater London. I've mapped these points onto a map of London's waterways and bodies and then used the distance matrix function to measure the distance of each mistletoe point to its nearest waterways vector node, and then calculated the mean distance from water.

I have then generated a set of random points across the same area and repeated the same process giving me a mean distance from water for the random points.

I'd like to repeat the process generating 10,000 separate sets of random points and perform a distance matrix for each set, and calculate its mean distance from water.

Is it possible to automate this in QGIS?


texts: Is co kriging appropriate for interpolation of three spatial co variates?

texts:
I want to estimate an interpolation model(s) for daily minimum temperature, daily maximum temperature and daily rainfall. I have samples of approximately the same size for all three variables at different locations although for certain days and locations only 1 or 2 of the 3 covariates are sampled . All three variables have some correlation, especially minimum temperature and maximum temperature. I read that co-kriging is especially adequate in cases where:

one wishes to interpolate one variable that is relatively sparsely sampled
one has available another covariate which is more densely sampled

This is not my case. Rather I have available samples for three covariates and want to interpolate all three. I have the following questions:

Is it better to develop a separate ordinary/universal kriging model for each variable or to develop a co-kriging model for all three.
I am using days and locations with no missing values to estimate the parameters of the interpolation model , but want to use it to interpolate missing values in different days and locations. Is it possible to do this in cases where:

certain locations have all three variables missing
certain locations have one o two out of the three variables missing

I am working in R in case there are any special considerations to heed.


texts: I have got myself a dxf which is built with units in the wrong scale. How do I upscale that in QGIS?   Note it does open in QGIS, just in the wrong map units and location.


texts: WMS GetFeatureInfo for multiple overlapping layers at the same time (Geoserver)

texts:
I have two layers served via GeoServer: "locations and "forests".

I'm using GetFeatureInfo request to query the layers like this:

mysite.net/geoserver/project/wms?REQUEST=GetFeatureInfo&SERVICE=WMS&SRS=EPSG%3A4326&STYLES=&TRANSPARENT=true&VERSION=1.1.1&FORMAT=image%2Fpng&BBOX=16.52648448944092%2C45.730599026324924%2C16.569614410400394%2C45.7570142556995&HEIGHT=882&WIDTH=1005&LAYERS=project%3Alocations%2project%3Aforests&QUERY_LAYERS=project%3Alocations%2Cproject%3Aforests&INFO_FORMAT=text%2Fhtml&X=383&Y=390

This works fine, but in cases where layers overlap, I get only the response for first layer (locations), and I need featureInfo for both layers. What am I missing?


texts: Structure From Motion - A priori estimation of picture requirements

texts:
Does anyone know how, or even if it is possible, to estimate the necessary number (and placement) of photos to achieve a certain final digital elevation model resolution (X points per square meter) generated using structure form motion?

I would like to do yearly follow-up surveys of a Brook that we have a longitudinal profile of, and have not been able (so far) to find any literature that discuss optimal picture number, placement, resolution, overlap, etc., given  a certain terrain X with characteristics Y, over a length Z.

Any ideas?


texts: ArcMap Crashes because of increasing memory usage

texts:
I have this code and I tried to delete all the layers it creates to enhance it's performance... However ArcMap crashes after about 200+ files.

I am new to the subject but I checked the internet for solutions and none of them worked (e.g. cleared cache). I would be really grateful if I get any help. Please note that lyr[number] depends on the rank of the layer in the TOC.

# Import arcpy module
import arcpy
arcpy.env.workspace = "E:\\Directory"
arcpy.env.overwriteOutput = True
fcs = arcpy.ListRasters()
for fc in fcs:
    # Local variables:
    ascii_input= str("E:\\Directory\\"+fc)
    colors = "colors"
    ascii_out = "ascii_out"

    # Process: Build Pyramids And Statistics
    x1=arcpy.BuildPyramidsandStatistics_management(ascii_input, "INCLUDE_SUBDIRECTORIES", "BUILD_PYRAMIDS", "CALCULATE_STATISTICS", "NONE")

    # Process: Define Projection
    x2=arcpy.DefineProjection_management(ascii_input, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
    arcpy.ClearWorkspaceCache_management()
    # Process: Make Raster Layer
    arcpy.MakeRasterLayer_management(ascii_input, ascii_out, "", "-97.154167 40.133331 -90.207778 44.53785", "")
    del(x1)
    del(x2)
    # Process: Apply Symbology From Layer
    arcpy.ApplySymbologyFromLayer_management(ascii_out, colors)
    ############### ##################################################################

    mxd = arcpy.mapping.MapDocument("CURRENT")
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        x=str(fc)
        if elm.text == "344":
            elm.text = x[2:4]+"/"+x[4:6]+"/20"+x[0:2]+" "+x[6:8]+":"+x[8:10] #change text to the date time
    arcpy.RefreshActiveView()
    lyr=arcpy.mapping.ListLayers(mxd)
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    arcpy.mapping.RemoveLayer(df,lyr[8]) #make sure it's the correct layer
    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToJPEG(mxd, r"E:\\IFC4GISJPEG\\"+str(fc)+".jpg")
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if elm.text == x[2:4]+"/"+x[4:6]+"/20"+x[0:2]+" "+x[6:8]+":"+x[8:10]:           #change text from date time back to 344 (a random value I chose)
            elm.text = "344"
    arcpy.RefreshActiveView()
    lyr=arcpy.mapping.ListLayers(mxd)
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    arcpy.mapping.RemoveLayer(df,lyr[7]) #make sure it's the correct
    arcpy.ClearWorkspaceCache_management()


texts: How to get the command line of QGIS algorithm "Merge vector layers"?

texts: To insert this function to my program I want to know what is command line executed by the algorithm "Merge vector layers" using QGIS GUI.


texts: How to connect QGIS (linux) to SQL Server

texts:
I tried every thing I can think of, official Microsoft ODBC Driver for SQL Server, FreeTDS, etc.... I'm able to connect to the SQL Server DB from the terminal (isql,sqsh) I can use PostgreSQL's Foreign Data Wrapper to connect PG and MSSQL, but I can't connect QGIS. I get the same error:

Driver not loaded Driver not loaded

How do I connect linux QGIS to SQL Server?


texts:
I am trying to build a qgis-plugin for a raster format (say .rst extension), I am able to build it on windows but facing issues on Linux.

I am running "make plugin" command from frtms/rst folder, it runs fine but don't creates a gdal_RST.so file. Prior to this I already ran configure and make command.

Am I missing anything. This is how plugin tag looks like in my GNUmakefile:

$(PLUGIN_SO):   $(OBJ)

        $(LD_SHARED) $(LNK_FLAGS) $(OBJ) $(CONFIG_LIBS) $(EXTRA_LIBS) \
            -o $(PLUGIN_SO)


texts: Mapping attributes from input shapefile to output CAD file using FME?

texts:
I'm trying to do a simple data conversion using FME where I take a point shapefile and convert it to CAD (DWG).

Even though I have connected the attributes I want mapped into the output file, the resulting CAD files ends up with no attributes. If I use 'Extended Entity Data' as my Attribute Output parameter I get several attributes but not these 3 three user attributes. I'm new to FME. My input data has a Z value attribute but is not actually Z- Aware.


texts: Calculating geometry size as a circle, and combining that data

texts:
Given a way geometry (from OSM), I would like to convert them to "circles" - a center point + radius guaranteed to contain the whole object. Plus I would need a set of functions to perform various union operations:

calculate centroid (already done using shapely)
calculate smallest radius that encloses the whole way
be able to combine two circles, produce a new center+radius circle that contains both.
be able to combine one circle with one point (find a new center + radius that contain previous circle and the node)
(simple case) combine two nodes into a circle

Is there a (preferably Python) library to do this?  Ideally the library should calculate surface distance rather than direct distance for the radius.


texts: Using COSI Corr Interior Orietation Fiducial Marks in UAV Imagery?

texts:
I am working on COSI Corr an ENVI Plugin, i am using it for Landslide Displacement using UAV Imageries.

How do I select fiducial marks?

I dont know what does it mean by "X loc (pix) and Y loc (Pix) ".


texts: Spatial Analysis with Factors/Classifications

texts:
When doing analysis on spatial data with associated information that is numeric, it is possible to use techniques such as Kriging or Interpolation or Cluster Analysis to find out interesting relationships in your data. Then you can visualise using your choice of poison/GIS. Okay.

Now, I want to do the same type of things with spatial data that has associated non-numeric data, i.e. a "Factor" in R or a "Classification". However I cannot find techniques that allow one to do this.

Does anyone know of methods that can be used to perform this type of analysis? It would also be a bonus you know whether they can be performed in R or QGIS.

I would really love a model to do interpolation with a non-numeric attribute. The set I would use this interpolation on would be lithology or other geological data.


texts: Reducing hourly rainfall data to monthly data for particular period of time and making a timeseries chart

texts:
I am working with GSMaP image collection (hourly rainfall data) in Google Earth Engine. I would like to make a script which is calculating and making a chart of sum rainfall in particular location during the March-May over the last 10 years. I made a script based on my skills in the platform, please find it below. The made script runs very slow and at the end gives error: "Error generating chart: Computation timed out". How can I optimize the script and solve this issue?

var point = ee.Geometry.Point([68.27212640755545, 40.007493398363955]);

var img_1 = ee.ImageCollection("JAXA/GPM_L3/GSMaP/v6/reanalysis")
    .filterDate('2000-03-01', '2014-03-01').select('hourlyPrecipRateGC');

var img_2 = ee.ImageCollection("JAXA/GPM_L3/GSMaP/v6/operational")
    .filterDate('2014-03-01', '2030-12-31').select('hourlyPrecipRateGC');

var imgcol_2008 = img_1.merge(img_2)
    .filterDate('2008-03-01', '2008-05-31').sum();

var imgcol_2009 = img_1.merge(img_2)
    .filterDate('2009-03-01', '2009-05-31').sum();

var imgcol_2010 = img_1.merge(img_2)
    .filterDate('2010-03-01', '2010-05-31').sum();

var imgcol_2011 = img_1.merge(img_2)
    .filterDate('2011-03-01', '2011-05-31').sum();

var imgcol_2012 = img_1.merge(img_2)
    .filterDate('2012-03-01', '2012-05-31').sum();

var imgcol_2013 = img_1.merge(img_2)
    .filterDate('2013-03-01', '2013-05-31').sum();

var imgcol_2014 = img_1.merge(img_2)
    .filterDate('2014-03-01', '2014-05-31').sum();

var imgcol_2015 = img_1.merge(img_2)
    .filterDate('2015-03-01', '2015-05-31').sum();

var imgcol_2016 = img_1.merge(img_2)
    .filterDate('2016-03-01', '2016-05-31').sum();

var imgcol_2017 = img_1.merge(img_2)
    .filterDate('2017-03-01', '2017-05-31').sum();

var imgcol_2018 = img_1.merge(img_2)
    .filterDate('2018-03-01', '2018-05-31').sum();

var imgcol_2019 = img_1.merge(img_2)
    .filterDate('2019-03-01', '2019-05-31').sum();

var imgcol_2020 = img_1.merge(img_2)
    .filterDate('2020-03-01', '2020-05-31').sum();

var image = ee.ImageCollection([imgcol_2008, imgcol_2009, imgcol_2010, imgcol_2011, imgcol_2012,
        imgcol_2013, imgcol_2014, imgcol_2015, imgcol_2016, imgcol_2017, imgcol_2018,
            imgcol_2019, imgcol_2020]);

var chart = ui.Chart.image.series(image, point, ee.Reducer.mean(), 1000);
  chart.setOptions({
    title: 'March-May Rainfall over years',
    vAxis: {title: 'Rainfall (mm)'},
    hAxis: {title: 'date', format: 'yy'},
  });

  print(image)
  print(chart)


texts: Accessing each raster in stack with loop using R

texts:
I have a stack of rasters and I want to access each layer, combine each later with another file and finally export the resulting file as a csv file.
Is there any way to do this with a loop?
So far I have done this but I do not get a file for each raster layer:

for (i in 1:nlayers(sms))
  {
  x<-sms[[i]]
  assign(x,raster(x))
  y<-extract(x,zipcoords)
  write.table(y,file=paste0("C:home",i,".csv"),sep=",",row.names=F)
  }


texts: Adding address block on Google Map Maker

texts:
If it is possible to mark a block of addresses for example B-block, xyz colony, which is a logical block of addresses on Google Mapmaker, how do I do it?

I want to add a polygon feature and select the block but in the options there is no 'block' choice or anything similar.


texts: Where can I find information about the development history of QGIS?

texts:
I know very little about early versions of QGIS (far far away from now..) and i'm feeling curious :-)

I am looking for various information about the past of QGIS, including:

Do you know in what year the software was created?
By whom? (Was it a single person or a group? Do you have some names? Nationality?)


texts: Adding a layer can sometimes load raster footprint!

texts:
Just trying to refine my ArcObjects code to add a layer into ArcMap based on passing the fully quailified (String) name of the layer - problem is, if it's a raster, it will load in the vector footprint of the raster which is no good.  I want to load the raster in, not it's footprint.

I can't see anything obvious that allows you to test for rasterness, other than swapping the test in my code around - i.e testing for raster first then vector.  Is there anything I am missing, is there a more elegant code out there to achieve this?

CODE

public static ILayer GetLayerFromSourceName(String name)
        {
            ILayer esriLayer = null;
            IWorkspace2 wsp2 = (IWorkspace2)Workspace;

        try
        {
            if (wsp2.get_NameExists(esriDatasetType.esriDTFeatureClass, name))
            {

                IFeatureLayer featureLayer = new FeatureLayerClass();
                featureLayer.FeatureClass = GeodatabaseManager.FeatureWorkspace.OpenFeatureClass(name);
                esriLayer = (ILayer)featureLayer;
            }
            else if (wsp2.get_NameExists(esriDatasetType.esriDTRasterDataset, name) ||
                wsp2.get_NameExists(esriDatasetType.esriDTRasterCatalog, name))
            {
                //This works only for remote raster datasests.
                IRasterDataset rasterDataset = null;
                IRasterLayer rasterLayer = new RasterLayerClass();
                rasterDataset = RasterWorkapaceEx.OpenRasterDataset(name);
                rasterLayer.CreateFromDataset(rasterDataset);

                esriLayer = (ILayer)rasterLayer;
            }
        }
        catch (Exception ex)
        {
            Debug.WriteLine(ex.Message);
        }
        return esriLayer;
    }


texts:
My organization has a web application for displaying global and regional climate data and we're considering improvements to it.  Our leading idea is to create an OpenLayers-based application to display climate rasters over base layers like OSM or Google Maps.  The map may also include some other application-specific vector layers and a color bar or other legend features describing the data.  Near as I can tell, this all seems do-able with OpenLayers.

One requirement for the project, though, is to be able to export (ideally from the web-app) publication quality (more-or-less) maps.  A typical use case would be something like this:

A user, using the base layer as a guide, navigates to their area of interest
The user selects a climate parameter (e.g. mean temp for the 2040s) or set of parameters
The user hits the "export" button and then downloads a good quality map which includes all of the elements described in the first paragraph.

I haven't been able to find anything in OpenLayers about rendering to an image.  Does anyone know whether this is possible?  And if not what other approaches for rendering could we take?  I would suppose that we would need to pass all of the map parameters (bbox, projection, included layers, etc.) out to some external rendering engine; what open source rendering engines are out there?


texts: Tool script multi-value parameter, Input raster does not exist

texts:
How do I adjust my script to accept a raster with a filename like the example below when used as a tool in ArcGIS with MulitValue set to Yes?

  "008. STUDY AREA MAP_rectified_NAD83_Z12.tif"

With files named in the same way as above the script works in python with two rasters separated with a semi-colon(;).

It also works as a tool when the input parameter (Raster Dataset) has MultiValue set to no.

If I use input file with more conventional names (just numbers and letters) the MultiValue option works correctly.

However, with MultiValue set to yes it fails at the first GetRasterProperties_management task with:

  ERROR 000865: [filepath\filename] does not exist

Unfortunately it is not possible to rename the files. Script below;

ScannedMaps = r'C:\ArcGIS\GeoRefCentrepoint\Data\010. STUDY AREA MAP1_rectified_NAD83_Z12.tif;C:\ArcGIS\GeoRefCentrepoint\Data\009. STUDY AREA MAP2_rectified_NAD83_Z12.tif'

        for image in ScannedMaps.split(;):
            print 'Starting: ',image
            arcpy.AddMessage('Starting: '+image)
            #Capture extent of georeferenced scanned map
            MaxX = float(arcpy.GetRasterProperties_management(image,"RIGHT").getOutput(0))
            MinX = float(arcpy.GetRasterProperties_management(image,"LEFT").getOutput(0))
            MaxY = float(arcpy.GetRasterProperties_management(image,"TOP").getOutput(0))
            MinY = float(arcpy.GetRasterProperties_management(image,"BOTTOM").getOutput(0))
            print '- Extent captured'
            arcpy.AddMessage('- Extent captured')

            #Calculate centre point
            centreX = (MaxX + MinX)/2
            centreY = (MaxY + MinY)/2
            print '- Centre point calculated'
            arcpy.AddMessage('- Centre point calculated')

            #Add centre point to feature, change file path, populate fields
            NewFilePath = image.replace('S:','\\\UNCNAME\FOLDER\SUBFOLDER')
            FileName = os.path.basename(image).split('.')[0]
            NewRow = [(FileName, NewFilePath, KeyWords, (centreX,centreY))]
            cursor = arcpy.da.InsertCursor(Feature,['DOCUMENT_NAME','URL','KEY_WORDS','SHAPE@XY'])
            print '- Attributes created'
            arcpy.AddMessage('- Attributes created')
            for data in NewRow:
                cursor.insertRow(data)
                print '- File:',FileName,'added successfully\n'
                arcpy.AddMessage('- File '+FileName+' added successfully')
            del cursor
    print 'Complete'

Update -

Output of when the script runs successfully from in PythonWin, included print ScannedMaps statement as requested:

Starting tool:

ScannedMaps = C:\ArcGIS\GeoRefCentrepoint\Data\010. MAP_rectified_NAD83_Z12.tif;C:\ArcGIS\GeoRefCentrepoint\Data\009.  MAP_rectified_NAD83_Z12.tif
Starting:  C:\ArcGIS\GeoRefCentrepoint\Data\010. MAP_rectified_NAD83_Z12.tif
- Extent captured
- Centre point calculated
- Attributes created
- File: 010 added successfully

ScannedMaps = C:\ArcGIS\GeoRefCentrepoint\Data\010. MAP_rectified_NAD83_Z12.tif;C:\ArcGIS\GeoRefCentrepoint\Data\009. MAP_rectified_NAD83_Z12.tif
Starting:  C:\ArcGIS\GeoRefCentrepoint\Data\009. MAP_rectified_NAD83_Z12.tif
- Extent captured
- Centre point calculated
- Attributes created
- File: 009 added successfully

Complete


texts: Classify forest with NDVI from Landsat 8 and 7

texts:
I'd like to calculate the NDVI for multi-date satellite images of Landsat (7,8) to perform a change detection on the deforestation of the Amazon rainforest. I already acquired surface reflectance products of Landsat 8 and 7. After calculating the NDVI, I'd like to reclassify the NDVI raster with the help of a threshold into at least forest and non-forest classes.
I tried to calculate the NDVI for two images (2000 and 2015) and found that I couldn't use the same threshold value to identify forest from non-forest in both images. So, my question is:

1) Does the threshold need to be the same for both images ? (I assume yes, but I'm a beginner.)


texts: multispectral images (Band Sequential format)

texts: How can I create a multispectral image file (in band sequential (BSQ) format), from several 1-band images, using GDAL libraries? I work with C++ language. Thanks for your attention.


texts: Calculate net elevation gain from a vector?

texts: I have a vector data (path/route) with points equally spaced with elevation data (in the attributes table). The elevation data was extracted from a raster. I would like to calculate net elevation gain from the starting point to the ending point. Is there a way to do this in QGIS?


texts: Comparing field of a table to calculate another field using ArcPy cursor?

texts:
How can I compare one field of a table to calculate another field? I have a table containing 'SL', 'Name' and 'Value' fields. I need to calculate the total value each name have. Example is shown in the picture below.
It will be great, if it would be Python script.


texts: Adding Values together using FindLabel (VBscript) function of ArcMap Label Expression?

texts:
I would like to add values from multiple fields in a shape file and display the result as a label

This is what I'm doing currently::

if ( [SwineGesBo] + [SwineSowLt] > 0) then
    output = output & [SwineGesBo] + [SwineSowLt] & "Sows"

Unfortunately the code here only puts the two value together and does not sum them.

How can I do this?


texts: Unable to load GeoJSON from external file in Internet Explorer 11

texts:
I created a web map in Leaflet using QGIS and the QGIS2Web Plugin. I want to be able to copy and paste this file onto a tablet, and also keep this file on the local network drive for my coworkers to access. I noticed that QGIS allows me to save a file as GeoJSON and the only difference between that and the JS files in the data folder is the extra code var variable =at the beginning of the latter.

I tried saving the GeoJSON into the data folder and loading it using an AJAX request and while it worked fine in Firefox as well as Internet Explorer 11 when I went to the F12 Developer Tools and put it in IE10 mode, it didn't load in IE11. The F12 developer tools told me that there were extra commas, but after deleting those it told me it didn't support the measure and search tools (which is weird as they work fine in the original html file), and then after I commented those out it would crash and give me a "this file stopped loading due to a long-running script" message.

Ultimately all I want to do is save myself the extra steps of saving my GeoJson file, opening it in Notepad, adding var variable = to the beginning of the code, and saving it as a JS file. Is there any way that I can do this without using an AJAX request, as this doesn't seem to work in Internet Explorer 11? I feel like this should be an easy thing to do.


texts: Using Symbols from Feature Layer for Graphics Layer in ArcGIS API for JavaScript?

texts:
I have a feature layer being drawn from a map service. I have a graphics layer that is the result of a query (happens to be on the same feature layer). I wish to apply the same symbols from the Feature Layer to the new graphics layer produced by the query but have had no such luck.  Relevant code bellow:

var featureLayer = new FeatureLayer(http://........);

function handleQuery(queryRslt) {
    var queryGraphicsArray = queryRslt.features;
    var featureLayerSymbol = featureLayer.renderer.getSymbol();
    for (var i = 0 ; i < queryGraphicsArray.length ; i ++) {
        queryGraphicsArray[i].symbol = featureLayerSymbol;
        map.graphics.add(queryGraphicsArray[i]);
}

I can confirm that the query is in fact functioning properly. I can also confirm that the symbol appears to being changed using console.log and debugging tool. However nothing is added to the map.

The symbol is a picture symbol added inside of arcMap before publishing to AGS. It is a ping image.

Any thoughts as to where i am going wrong?


texts: Making polygon boundaries of different feature classes coincident?

texts:
I have two polygon feature classes:

Cadastral layer and
another polygon layer "zones".

In both feature classes the polygon boundaries are coincident inside their feature but not between the features.

What I want to do is leave the cadastral layer boundaries fixed and snap the 'zones' layer boundaries to them as close as possible within a distance of 2 metres to remove all the small silvers.

I have tried snapping zones to cadastral with a distance of 2 metres and that worked for some boundaries but not all slivers 2m of less. Why do some boundary lines snap to each other when inside the 2m distance and some do not?

I have also tried integrate where you can rank one class over another and generally it works well for distances 1 metre or less but over 1m can even move the rank 1 class boundaries when snapping the boundaries together which is not suitable.

What tool or combination of tools can I use to snap the boundaries of the different features together when they are within 2 metres of each other and boundaries outside of 2 metres stay as they are?

The only strong rule I want to use is that the cadastral layer boundaries cannot move at all.


texts: Module Catchment area (Parallel) missing from QGIS processing?

texts: I have qgis 2.4.0 chugiak and have installed processing but in the box of tools not exist the module catchment area (parallel), how to install this module?


texts: OpenLayers ArcGIS93Rest layer not showing up

texts:
So this seems like it should be super obvious, but I can't for the life of me figure out why I can't get any ESRI layers to show up on my map. I've tried different layers and messing around with the options, but all that happens is that both layer names show up in the layers panel, but only the Google layer actually shows up on the map. Here is the code I'm using:

function addOpenLayers() {
    var mapcenter = new OpenLayers.LonLat(-90.8, 46.45)
    var base = new OpenLayers.Layer.Google("Google Satellite", {type: google.maps.MapTypeId.SATELLITE, numZoomLevels:22});

    var olOptions = {
        controls: [new OpenLayers.Control.PanZoomBar, new OpenLayers.Control.LayerSwitcher, new OpenLayers.Control.Navigation],
    }
    var olMap = new OpenLayers.Map("map", olOptions);

    olMap.addLayer(base);

    olMap.setCenter(mapcenter.transform(
        new OpenLayers.Projection("EPSG:4326"),
        olMap.getProjectionObject()
    ), 10);

    //This is the part that doesn't work:
    var esriRoadsUrl = "http://server.arcgisonline.com/ArcGIS/rest/services/Reference/World_Transportation/MapServer/export";
    var esriRoads = new OpenLayers.Layer.ArcGIS93Rest("Roads", esriRoadsUrl, {layers: "0"});
    esriRoads.isBaseLayer = false;

    olMap.addLayer(esriRoads);
};


texts: I have created the snow cover map using the NDSI technique. I want to create a graph showing the snow's areal extent with progressing time. How can I create it?


texts: Is it possible to limit export feature class in ArcSDE to only some users?

texts: I'm trying to manage user and their privileges in ArcSDE and sql server 2012. I can privilege the select and edit a layer for users. I can use just two privileges but i want to permit export feature class capability for some users.Is it possible to permit and revoke export feature class for users ?


texts: How to re-grid a raster using area-overlap weighted averaging?

texts:
I have two rasters r and r2 with different resolutions. I want to re-sample (re-grid) r2 to be the same as r using area-overlap weighted averaging:If an r grid cell is covered by more than 1 r2gridcell, then the value assigned is the added value of each of those r2 cells weighed (multiplied) by the % of the r cell they cover.

There are some available methods in gdal but I think none of them matches this method. Available methods are:

           near, bilinear, bilinear resampling,cubic,cubic resampling,cubicspline, average ...etc

Has anyone done this before?


texts: Finding parcels which combine to create largest area

texts:
I have a shapefile of 500 polygons and I need to find the 22 neighboring polygons that when combined will produce the largest continuous area.

Does this involve the Dissolve tool?

What do I dissolve the polygons by in order to get the largest area?


texts:
I made a silly mistake: when creating linestrings in my PostGis database i inverted the lon and lat values. Which isn't a problem for postgis itself, but for my visualization and intersections. ;)

Of course i only recorgnized this after importing lots of stuff. Instead of simply re-running the import i'm curios if there is a way i could swap the X/Y-Coordinates of every point of a linestring?

e.g.:

LINESTRING(lat1 lon1, lat2 lon2, lat3 lon3); // what i got
SomeMagicMethod();
LINESTRING(lon1 lat1, lon2 lat2, lon3 lat3); // what i want


texts: Import shapefiles into PostGIS using Python with psycopg2

texts:
I want to import shapefiles to PostGIS using Python with arcpy & psycopg2 to connect to PostGIS.
I don't want to use the existing shapefile loader nor ogr, just write my own import tool.

My steps and where I got stuck:

Define shapefile
Select column names from shapefile into a list
Create SQL CREATE TABLE query with column names of shapefile
Insert values from shapefile into PostGIS database

I got stuck at point 4 since I don't really know how to extract the values from the shapefile and insert it into the database taking account of the correct column.

Any idea?

Code example:

import arcpy
import psycopg2

arcpy.env.workspace = r"path\to\data"

shp = "BRD.shp"
shp2 = shp[:3]

desc_shp = arcpy.Describe(shp)

#set variables and create list wo OID and Shape attributes
if desc_shp.hasOID:
    print desc_shp.OIDFieldName
    shp_OID = desc_shp.OIDFieldName

print desc_shp.shapeFieldName
shp_geom = desc_shp.shapeFieldName

fields = desc_shp.fields
field_names = []
for field in fields:
    if not (field.name == shp_OID or field.name == shp_geom):
        field_names.append(field.name)
new_list = field_names
print new_list

sc_da_cursor = arcpy.da.SearchCursor(shp,shp_geom)
for sc_da_row in sc_da_cursor:
    print sc_da_row[0]

# PostGIS part after here

user = "postgres"
pw = "postgres"
db = "pypg"
host = "localhost"
port = 5432

conn_obj = psycopg2.connect(database=db, user=user, password=pw, host=host, port=port)
cur_obj = conn_obj.cursor()

sql_createtable = "CREATE TABLE " + shp2 + "(" + shp_OID + " serial PRIMARY KEY," +
shp_geom + " geometry," + new_list[0] + " VARCHAR(100)," +  new_list[1] +
" VARCHAR(100))"

sql_insert = "INSERT INTO " + shp2 + "(" + shp_geom + ") VALUES(" + str(sc_da_row[0]) +
")"
cur_obj.execute(sql_insert)

conn_obj.commit()
conn_obj.close()
cur_obj.close()


texts: Publishing feature service from local file geodatabase using ArcGIS for Server?

texts:
How is it possible to publish the feature service from a local file geodatabase i.e. not using an Enterprise geodatabase such as Oracle via ArcSDE.

When i try to create a feature service it gives me error:

  Feature Service require a registered database

and i have no database installed. i am copying my geodatabase to ArcGIS Server.

When i try to place my file geodatabase on network that is shared location it gives me error that:

  Feature Data Source is not supportted

Software that I am using:

ArcGIS for Desktop 10.4.1
ArcGIS for Server 10.4.1
Windows 7 64 bit


texts:
I have a code that calculates NDVI for an imageCollection and then is supposed to display those rasters. The rasters are filtered by their size so only rasters that have more than a certain amount of pixels are being added to a new imageCollection and then displayed.
The problem is that the images cannot be displayed and I get the next error message:

  Layer error: Element.get, argument 'property': Invalid type. Expected:
  String. Actual: Integer

The graphs are displayed but the images are not and I don't understand why would the images need to be string, when did they changed and how to change them back to integer.

This is my code:

 * Function to mask clouds based on the pixel_qa band of Landsat SR data.
 * @param {ee.Image} image Input Landsat SR image
 * @return {ee.Image} Cloudmasked Landsat image
 */
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2).divide(10000)
  .copyProperties(image, ['system:time_start']);
};

//Create LANDSAT7 dataset

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filterDate('1999-01-01', '2013-04-29')
                  .select('B1','B2','B3','B4','pixel_qa')
                  .filterBounds(geometry)
                  .map(cloudMaskL457);

//RGB visualization

var visParams = {
  bands: ['B3', 'B2', 'B1'],
  min: 0,
  max: 1,
  gamma: 1.4,
};

//clip the dataset according to the geometry
var clippedCol=dataset.map(function(im){
  return im.clip(geometry);
});

// Get the number of images.
var count = dataset.size();
print('Count: ',count);
//print(clippedCol);
print(dataset,'dataset');

//function to calculate NDWI in LANDSAT7
var addNDVI = function(image) {
  var NDVI = image.normalizedDifference(['B4', 'B3'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(NDVI);

};

//NDWI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var NDVIcolor = {
  min: -1,
  max:1,
  palette: ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'],

};

var ndviWithCount = withNDVI.map(function(image){
  var countpixels = ee.Number(image.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: geometry,
  crs: 'EPSG:4326',
  scale: 30,
  }).get('NDVI'));

  return image.set('count', countpixels);
});

print(ndviWithCount, 'ndviWithCount');

//filter between a range
var filter = ndviWithCount.filter(ee.Filter.rangeContains(
          'count', 220000, 446224))
print(filter, 'filtered');

var max = ndviWithCount.reduceColumns(ee.Reducer.max(),  ["count"])
print(max.get('max'));

var OrderList=[0,1,2,3,4,5,6,7,8,9,10]

for (var i in OrderList) {
  var image = ee.Image(filter.get(OrderList[i]));
  var toexport=image.visualize(NDVIcolor).addBands(image);

  // do what ever you need with image
  Map.addLayer(image, NDVIcolor, i);

 }

print(ui.Chart.image.series(filter, geometry, ee.Reducer.mean(), 30));
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.max(), 30));
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.min(), 30));
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.stdDev(), 30));


texts: Implement Layers of Details (LODs) for custom basemaps using ArcGIS Javascript API 3.17

texts:
I am trying to implement layers of details (LODs) on a few custom basemaps using the ArcGIS Javascript API 3.17. One is for less detail, and the other two used for more detail, but covering different zones. I'm having issues, and the page isn't loading however (for Firefox and IE). Can anyone see anything I'm not? Note, that I've used some sample services here. Additionally, I've added some dynamic feature layers in my real implementation, that I've not included here.

var map;
var extent;
var baseMapLayerPPT;
var baseMapLayerPoly;
var customLods = [];
var loadCount = 0;

require([
  "esri/geometry/Extent","esri/map", "esri/dijit/Search", "esri/layers/ArcGISDynamicMapServiceLayer","esri/layers/FeatureLayer", "esri/InfoTemplate", "dojo/domReady!"
], function (Extent, Map, Search, DynamicMapServiceLayer, FeatureLayer, InfoTemplate) {

  baseMapLayerPoly = new esri.layers.ArcGISTiledMapServiceLayer("https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer", {
    displayLevels: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
  });

  baseMapLayerPoly.on("load", addLods);

  baseMapLayerPPT = new ArcGISTiledMapServiceLayer("https://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer", {
    displayLevels: [11, 12, 13, 14, 15 ],
    //opacity : 0.75
  });

  extent = new esri.geometry.Extent({
    "xmin":-91070.37,"ymin":7988806.64,"xmax":301214.62,"ymax":8227343.05,
    "spatialReference":{"wkid":3297}
  });

  // popuplate an array with zoom levels
  function addLods(evt) {
    customLods = customLods.concat(evt.layer.tileInfo.lods);
    loadCount++;
    if (loadCount === 2) {
      initMap();
    }
  }

  // create the map and use the custom zoom levels
  function initMap() {
    map = new Map("mapDiv", {
      center: [-150, -17],
      zoom: 9,
      lods : customLods
    });
    map.on("extent-change", changeScale);
    map.addLayers([baseMapLayerPoly,baseMapLayerPPT);
  }

  //Report which layer is being shown
  function changeScale(evt) {
    dom.byId("scale").innerHTML = "Level: " + evt.lod.level;
    if (evt.lod.level < 11) {
      dom.byId("visibleLayer").innerHTML = "Layer: Imagery";
    } else if (evt.lod.level == 11) {
      // both layers are loaded
      dom.byId("visibleLayer").innerHTML = "Both layers currently visible";
    } else {
      dom.byId("visibleLayer").innerHTML = "Layer: Streets";
    }
  }

});


texts: Inverse of Euclidean distance in ArcGIS Desktop?

texts:
Using ArcGIS Desktop I am trying to do the reverse / inverse of a euclidean distance. For example, I am doing distance to roads analysis. The output from the distance raster shows roads as 0 and areas closest to roads as 0.1, 0.2, etc and the areas furthest from roads as a very high value. However, I would like the raster to be the reverse of this - the areas closest to the roads should have the highest values and the areas furthest from the roads should have the lowest values.

How do I create a raster that reverses the number?

I would use reclass but this would be difficult to use as the values are continuous and not integers.


texts:
In Arc, I am used to setting a fixed scale for the Data Frame, so when you zoom in and out the objects get larger or smaller. Is there a way to set a fixed scale in QGIS? At the moment, if you zoom in and out, the lines and text stay the same size.
I normally use ArcGIS, but we only have a single licence, so we are doing some research to using QGIS as a standby for some of our work.


texts: How to generate non-merging buffers in QGIS?

texts:
I'm using QGIS to construct some fixed width buffers from a coastline shapefile. Is there a way to easily handle overlaps/ median lines? (I'm not sure if there's a particular term or phrase I should be searching for here.)

More detail:

I have a shapefile for several islands. I need to create a 12nm territorial sea coastal zone. This is straightforward with the buffer function.

However, if the distance between two islands is less than 24nm (12 + 12), the boundary should follow the median line.

Is there a way to calculate this in QGIS so that, instead of a merged polygon, I can generate two polygons (one around each island) with a boundary along the mid point.

For bonus points, I need to generate a second smaller coastal zone. This lies 4nm from the coast, or 1/3 of the distance to the territorial sea boundary described above. Is there a way to construct a buffer with a proportional distance?

I hope my explanations make sense. Please let me know if anything is unclear and I'll add more detail.


texts:
Trying to set up a defintion query from a connected database table in ArcGIS.

I am using this SQL (which runs properly in SQLDeveloper) but i keep on getting the error "The SQL expression has invalid syntax"

Any thoughts??

select distinct
DYN.%BR_REPORTING.GENERATION_TONS,
DYN.%BR_REPORTING.REPORT_CYCLE,
DYN.%BR_REPORTING.HANDLER_NAME,
DYN.%BR_REPORTING.HANDLER_ID
from DYN.%BR_REPORTING
where 1=1 and upper(DYN.%BR_REPORTING.REPORT_CYCLE) = 2017
order by
DYN.%BR_REPORTING.GENERATION_TONS,
DYN.%BR_REPORTING.REPORT_CYCLE,
DYN.%BR_REPORTING.HANDLER_NAME,
DYN.%BR_REPORTING.HANDLER_ID


texts: Select rows where column values don't match the domain

texts:
In ArcGIS Desktop, I can use the Validate Features option in the Editor Toolbar menu to validate feature class columns against their respective domains. More info here in the Validating Features section.

However, this functionality only works on spatial feature classes, not on non-spatial tables. Additionally, even if I were using a spatial feature-class, the feedback that the Validate Features function provides is too limited to be useful.

I want to select records in non-spatial tables where column values don't match their respective domains.

How can I do this?


texts: What is the advantage of using Postgres with Esri architecture?

texts:
One of my tasks is to transfer shapefiles and .gdb data into a PostgreSQL database - not a problem. However, the database is to be connected to ArcGIS Server so that the data can be used by an ArcGIS Online app. This requires the creation of a geodatabase within the PostgreSQL database so that ArcGIS can connect to it.

This seems inefficient: The workflow, as communicated to me, is to receive data (probably as a shapefile) > Insert into PostgreSQL > Use geodatabase within PostgreSQL to provide data to ArcGIS Server.

The only advantage I can think of is that using PostgreSQL is a way to allow for easier access by other software. However, if this database is going to be used by Esri architecture, why not use a geodatabase?


texts: Iterate Workspaces in Modelbuilder not executing all tools

texts:
I have a folder which contains a number of file gdbs. These gdbs all have the same structure i.e. a feature dataset FD_1 and several feature classes within the dataset (the two I am interested in are FC1 and FC2). The feature classes all have the same names within the gdbs.

I want to create tables from the feature classes, so my model does the following:
Iterate Workspaces produces Workspace (the path of the current gdb) and Name (the name of the current gdb). The Copy Rows tool is used to produce a table for FC1. The input table is specified using inline variable substitution

%Workspace%\FD_1\FC1

As is the output table

%Output%\table1_%Name%.dbf

The problem is, the Copy Rows process is not executing. The model simply iterates over the gdbs and then ends. When I check the Copy Rows window, it says that the input table does not exist. I tried setting Workspace as the Current Environment for the tool, but it didn't change anything.


texts: Are MapBox mobile SDK (iOS, Android) able to show raster tiles from binary images on memory?

texts:
I want to use MapBox mobile SDK (iOS, Android).
Actually I use them from Xamarin, but Xamarin SDK is just wrapper of Obj-C/Java SDKs, so I think there are no differences.

What I want to do is, showing raster tile maps from binary images on memory.
In my app, I have tile image data as binary in memory, like NSData or byte, so I show them on map as web Mercator raster tiles.

I can't find such interface on MapBox SDK's document.

Anyone know whether MapBox SDK can raster tiles from binary image on memory or not?


texts: Different ways to publish media-enhanced route/trail information to mobile users

texts:
I'm looking for some different approaches that a non-profit I'm helping could use to publish information about some tourism routes. Each route is 1-2 tracks, with some historical  background information in the form of text and pictures, plus some POIs representing specific historical details and logistics (eg, toilets, parking, cafes, hotels...) We'd like to be able to give this to people in a form they can use on mobile phones, while enjoying the trails - preferably available offline. Ideally, they would be able to see at a glance "20m on my left is a former rail station [see photo], and the nearest food is 3km further ahead up a steep hill". The view would (ideally) be linear, rather than just a top down 2D map view.

So, some approaches we've thought of so far:

Develop an app for each route
Find an app that does most of what we want, and provide content for it (everytrail.com is one example - would like to hear of others)
Provide track and POI information in KMZ form, choose a recommended app or two (eg, Locus), and give instructions on how to combine them.

We're also adding information to OpenStreetMap where appropriate.

So - I would love to hear of any other possible ways to tackle this. We may well pursue more than one simultaneously.

[Apologies if this isn't quite within the scope of GIS. It's not quite within the scope of SO, SU or Android either...]


texts: Using python QGIS network analysis library to calculate short path based on speed and distance?

texts:
I want to calculate and show short path from starting point to an end point in a qgis python plugin.
in my research i stopped by an article about QGis network analysis library; http://www.qgis.org/en/docs/pyqgis_developer_cookbook/network_analysis.html,
and i found this:
"To calculate edge properties programming pattern strategy is used. For now only QgsDistanceArcProperter strategy is available, that takes into account the length of the route. You can implement your own strategy that will use all necessary parameters. For example, RoadGraph plugin uses strategy that compute travel time using edge length and speed value from attributes."
I'm wondering if i can define my own strategy based on (speed&distance) meaning time cost using QgsDistanceArcProperter.

Can you explain to me if it is possible and how?


texts: Shortcut to copy current cursor coordinates in ArcMap

texts:
Is there a shortcut or an extension to ArcMap 10 that would allow me to quickly copy (to clipboard) the coordinates at the current map position/click?  I could implement such a functionality easily enough, however I'm hoping someone has already done so.

edit: Apparently the code is much simpler than I figured:

IPoint pt = ArcMap.Document.CurrentLocation;
Clipboard.SetData(DataFormats.StringFormat, pt.X + " " + pt.Y);


texts: View Cached Google Map in Open Data Kit's Collect Application

texts:
I am using Open Data Kit to collect data on an Android tablet.  The ODK form was built using the ODK Build application on-line.  The Build application allows the user to record a latitude and longitude and view the current GNSS location on a Google Map background-this works great if I have a wifi connection.  Since my tablet does not have a cellular connection I have cached the Google Maps tiles on to the tablet.  I can view the tiles when there is no wifi connection in the Google Maps app so I know the caching works.  I cannot view the cached map tiles while collecting data in the field while using the ODK Collect app without wifi.
  Since these Android devices store the last know location it might record/report an invalid location given GNSS availability. I would like to see the current location I am recording on a map to ensure that I am in the correct location.

Is there some way to see the cached map tiles in the ODK Collect app on the tablet?


texts: Is there a way to take control of the mouse pointer using python in QGIS?

texts: I would like to take control of the mouse pointer in QGIS (for example constraining it to a given extent...) Is there any way to do that using PyQGIS?


texts: Enable snapping options in QGIS-2.18 plugin using Python

texts:
I have a custom plugin to place point feature on another point feature.
So i want to enable/activate snapping options before i place that feature in my QGIS application using python.

So far i tried lot of answers that i found here

Set Snap Mode with QGIS

enable snapping in QGIS

but none works for me ...

Anyone help me please ..!!!


texts:
I have an asc type file which describes measurements of a gas across an area at 1km square resolution. It looks like this:

NCols 1725
NRows 2175
xllcorner -224131.189661
yllcorner 4892726.61578
cellsize 1000.0
NODATA_value -1.0
-0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0 etc

I can visualise this raster format in QGIS:

I now want to convert this to a vector form.

Ultimately I wish to convert this to a csv or json type file that contains the bounds for each 1km square section and the associated metric being visualised in the raster format.

I have tried to convert from raster to vector using gdal as an intermediate step in the right direction:

gdal_polygonize.py sample_file.asc

However, this resulted in:

"/Library/Frameworks/GDAL.framework/Versions/1.11/Programs/gdal_polygonize.py", line 186, in
    dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs )
    AttributeError: 'NoneType' object has no attribute 'CreateLayer'

Is there a more direct path for doing this?

Edit. Running gdalinfo file_name presents:

Driver: AAIGrid/Arc/Info ASCII Grid
Files: file.asc
Size is 1725, 2175
Coordinate System is `'
Origin = (-224131.189661000011256,7067726.615779999643564)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Corner Coordinates:
Upper Left  ( -224131.190, 7067726.616)
Lower Left  ( -224131.190, 4892726.616)
Upper Right ( 1500868.810, 7067726.616)
Lower Right ( 1500868.810, 4892726.616)
Center      (  638368.810, 5980226.616)
Band 1 Block=1725x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-1

From this I can calculate the bounds of each invidual pixel simply enough. However, the issue of taking the actual data field itself persists.


texts: Assigning names to features when exporting to KML in Google Earth Engine

texts:
I'm working in Google Earth Engine using Javascript api and I'm exporting my feature collection as a KML file to open in Google Earth Pro. One of the properties of each feature is its station_id. Right now, all the features are labelled "no name" when I open the KML in Google Earth Pro. I want the station_id to become the feature names when I open up the KML file in Google Earth Pro.
Is there a way to change the feature collection in GEE so that each feature has a name when I open it in Google Earth? Alternatively, is there a way to do this in Google Earth Pro?


texts: Building attribute table from raster in ArcMap?

texts:
I'm using ArcMap 10.1.

I have a NetCDF file containing the latitude, longitude and porosity for the entire US. I am trying to add another value (wilting point) which depends on the porosity.

I wanted to extract the netcdf to a table, add a new field and then use the field calculator to to fill the values on the new field. But when I import it, I only get one set of values on the table.

I tried approaching it the other way: extracting the NetCDF to raster (which works fine) and extracting the values from the raster using "build raster attribute table" but i get an error message saying "Only single band integer raster dataset is a valid input". it doesn't make sense as the raster is a single band raster (I double checked the properties> number of bands of the raster). HOw can I build a table containing all the entries in the NetCDF file?


texts: Difference between map layer and spatial dataset (e.g. feature class, shapefile)?

texts:
What is the difference between a feature class and a feature layer?

To me they seem very similar in concept.


texts: I'm new in using ESRI CityEngine 2011. I read it's tutorials. It was very interesting and helpful for but there is no references or tutorials for how to create CGA rules. How can I find any tutorials  about this subject?


texts: Why am I getting an error when using multiple logical filters for WFS?

texts:
I am using openlayers to query a geoserver WFS and am trying to use ol.format.filter.and to include multiple logical filters:

var featureRequest = new ol.format.WFS().writeGetFeature({
        srsName: 'EPSG:3857',
        featureNS: 'http://geoserver.sf.net',
        featurePrefix: 'sde',
        featureTypes: ['parcels_il'],
        outputFormat: 'application/json',
        filter: ol.format.filter.and(
            ol.format.filter.like('County', 'Henry*'),
            ol.format.filter.equalTo('Landcover Majority', 'Cultivated Cropland')

        )
      });

      // then post the request and add the received features to a layer
      fetch('http://localhost:8080/geoserver/wfs', {
        method: 'POST',
        body: new XMLSerializer().serializeToString(featureRequest)
      }).then(function(response) {
        return response.json();
      }).then(function(json) {
        var features = new ol.format.GeoJSON().readFeatures(json);
        vectorSource.addFeatures(features);
        map.getView().fit(vectorSource.getExtent());
        map.addLayer(vector);
      });

I can run each filter individually without a problem, but when I try to combine I get this error:

filter_test.html:1 Uncaught (in promise) SyntaxError: Unexpected token < in JSON at position 0

I've researched this but haven't figured it out yet. What is going wrong?

UPDATE: Upon further investigation it doesn't seem that the issue is with multiple logical operators in all cases, it seems to have issues with string fields where there are spaces in the values. I also get this error when trying to filter on numeric fields using lessThanorEqualTo and greaterThanorEqualto.


texts: Set minimum scale on QGis Print composer (atlas)

texts:
I have a template ready to be printed on Print composer tool of Qgis. I need to fix a minimum scale while figures are represented one by one (1:10000)
I would like to combine the scale tool controled by Atlas which sets a 10% margin of the figure, with a precondition of setting minimum scale at 1:10000.
Base cartography is 1:10000 and sometimes figure is really small so Atlas represent the figure with a margin of 10% so scale could be 1:2000 and resolution of base cartography is bad.


texts: QGIS OpenLayers OSM Plugin bypasses default SystemProxy

texts:
Tested with QGIS 2.8, 2.14 and 2.18 on same desktop Client.

QGIS Proxy-settings on network-settings are set to DefaultProxy.
OpenLayers plugin GoogleMaps works properly. Uses SystemProxy.

OpenLayers plugin OpenStreetMap allways tries to bypass the proxy and send requests directly into the internet. This is denied by firewall. WWW-access only allowed via proxy.

Tried to do diagnostics with Fiddler4. I see the requests for GoogleMaps passing Fiddler forwarded to SystemProxy. Everything fine.
QGIS-OSM requests do not appear in Fiddler. QGIS-OSM requests appear in Firewall log as denied requests.

So while QGIS-GoogleMaps works properly and QGIS-OSM does not use DefaultProxy, the question is: Is there any other setting in QGIS-Software where I can force OSM Plugin to direct the requests to Proxy?


texts: How to save many print view with the same settings?

texts:
I use QGIS 2.6.1 and I need to save more than 100 map shots (one for each Municipality) for future print and I need to save both the print settings and each shots (with positions and zoom levels).

With the print compositions I can save the settings but not the shots, savings the models I have the positions and the zooms but not the settings


texts: Creating admin roles with limited access to PostGIS

texts:
I'm trying to cram three different organisations into one server. Today they each have their own geoserver instance connected to separate PostGIS databases, which works fine. But now they also need to be able to access and manipulate their own data and I can't figure out how to set this up.

So in short: How do I set up three different users in PostGIS and grant them access to one specific database each?

Any suggestions welcome

/O


texts: Thiessen polygons of sub-districts to lie within district boundaries

texts: I'm using ArcGIS 10.1 and am fairly new to it. I have the latitude and longitude of sub-districts, and wish to create Thiessen polygons of these sub-districts. I also have shapefiles for districts (there are about 10 sub-districts within each district), and want the Thiessen polygon of each sub-district to lie completely within it's district, and not spillover to a neighboring district. Thanks!


texts: Difference between PostGIS ST_Intersects vs '=', QGIS and ArcGIS 'select by location'

texts:
I have two tables, 'protein' and 'data_aggreg', in a PostGIS database. I am simply trying to query the points that intersect (all of the points in 'data_aggreg' share exactly the same geometry as a portion of the points in 'protein' table).

This is typically a simple operation. However, this time PostGIS and QGIS return nothing and ArcGIS returns the correct selection of points. The data in both QGIS and ArcGIS are sourced directly from the same PostGIS database via their respective methods of connecting to external DBs (ie, no shapefiles, csv's, etc are being used for this analysis). Both tables are in projected coordinate system EPSG 26912 and both tables are single point geometry.

When I run ST_Intersects in PostGIS on the two geom columns nothing is returned. When I run the 'select by location' tool in QGIS nothing is returned (I've also tried spatial join with same effect).

Syntax for ST_Intersects, which returns nothing:

SELECT protein.protein, data_aggreg.pro_2016
from production.protein, scratch.data_aggreg
where ST_Intersects(geom26912, geom);

However when I run 'select by location' in ArcGIS I get the expected result. That is, the subset of points in the 'protein' table that share the same geometry with 'data_aggreg' are correctly selected. Also, when I change my PostGIS 'where' query from ST_Intersects to geom1 = geom2 I also get the expected result.

Syntax when using '=' in PostGIS, which returns the expected result:

SELECT protein.protein, data_aggreg.pro_2016
from production.protein, scratch.data_aggreg
where geom26912 = geom;

It is my understanding that if any part of features share any space then ST_Intersects will return 'true'. So it follows that if geom1 = geom2 then ST_Intersects should return 'true'. But this is not the case.

So why the difference in results within PostGIS when using ST_Intersects vs '=', and why the difference among PostGIS when using ST_Intersects, QGIS 'select by location', and ArcGIS 'select by location'?


texts: Is rid=2 significant in PostGIS raster structure?

texts:
It seems that in most, if not all, examples with raster appear to include the WHERE rid=2 statement in the query.

http://postgis.net/docs/RT_ST_Value.html

Is there something about the rast structure that requires this, or is this just an artifact of previous examples used to reduce the sample dataset?


texts: Creating space time cube gives ERROR 110015 Unable to read input?

texts:
Somehow the create space time cube doesn't work on my project on Philadelphia crime data. The input feature is the point feature "robbery", I named the output as robbery_STC.nc, the time field is one field in the attribute table named dispatched date, all other blanks are left default.
The error message is:

  ERROR 110015, Unable to read the input space time cube; and failed to
  execute (CreateTimeSpaceCube)


texts:
I'm new with GDAL ... and I'm trying to write some code where I have to read a .raw image and store it as an array of 32-bit floats. It seems that it is easy with python, I just have to use:

array = image.ReadAsArray().asType(float32)

But What I'm trying to know is if there is any way to have the same thing with the C API.


texts: Are there any document icons for the qgis mimetype?

texts:
In a default Ubuntu install (QGIS 1.8.0) no document icon is associated with the mimetype 'application/x-qgis-project'.  As a consequence every saved QGIS document is interpreted by the OS as being a text or xml document and the icon displayed gives no indication it has anything to do with QGIS.

I've googled around and can't find anything; I could make an icon set myself but wanted to see if I'm just missing something before doing it..


texts: What coordinate system/projection are these Fresno, CA data points?

texts:
I have some shapefiles which include various points and a corresponding street address, whose lat/lon I've manually found in Google Maps. There is no .prj file.

Google Maps Lat,Lon,       unknown coordinate system x, y, city, state
36.7448656,-119.7937943,   6329742.00029, 2155243.06469,   Fresno, CA
36.7341137,-119.8126966,   6324144.37016, 2151400.93634,   Fresno, CA


texts: Is it possible to learn ArcObjects without VB?

texts:
I have some text in arcmap active view and I need to Change the font face as the licence does not cover the font type. I am using python script for Automation but arcpy can not select the text elements from arcmap active view.

I was looking for solution and I found Arcobject can solve it. But I do not have Visual Basic. So, How can I learn ArcObject?


texts:
I have calculated a viewshed (ZVI) in OpenWind and opened the raster in QGIS

I have then converted the raster to a vector

When I open the attribute table there are 4682 rows (features) and there is a column titled 'DN' and each row contains a number between 1 and 10, also 4681 rows long

When I add an $area column to work out the area of the viewshed, each row says 'NULL'

I am simply trying to calculate the area of the viewshed in metres

I am not sure why there are so many rows or why the area returns as 'NULL'?


texts: Problem exporting registers with the land registry in DXF continuation (suite) in reprojection

texts:
I visualize (display) a municipality with the plugin register with the land registry in spatiality and in projection Lambert 93. I need the bottom on Autocad in projection CC47 thus I re-throw (re-plan) on QGIS by choosing as current SCR the CC47 (the reprojection in the air bustles). Up to there, everything works.

Then, in Project, I choose export DXF and I make my export. It sticks (there is a problem) when I visualize (display) the created file DXF because either I have no object inside, or I have all the municipality but in Lambert 93 (as I mark or not "exporter only the visible entities on the influence of the card (map)").

When I make "Project" then "Export DXF", I cannot choose the SCR ( CC47).


texts:
I have a list with lots of VRTs. Each VRT is the mosaic for a single band for a single band. What I want to do is to make a final TIFF with as many layers as bands. I'm trying to first make one VRT and then convert it to a TIFF in Python. But what I get is a TIFF with a single layer...

My list of VRTs (vrts_per_date) kind of looks like this:

[296B02.vrt', 296B10.vrt', 296B08.vrt', 296B12.vrt', 296B07.vrt', 296B03.vrt', 296B11.vrt'...]

And then I use:

single_vrt = gdal.BuildVRT("my.vrt" vrts_per_date)

And finally I want to make a TIFF with all bands, so I use:

final_tif = gdal.Translate("my.tiff", single_vrt)

But when I open it in QGIS it's a mosaic, but with just one band...
How can I get a TIFF as output with as many bands as VRTs in vrts_per_date?


texts: Buffer points by buffer amount in attribute value

texts:
I have a set of points in a shapefile with a real number attribute 'r' containing a radius (cp. screenshot below) and I would like to buffer these with amount of 'r'. As the result I expect a polygon shapefile with approximated circles of various radius, but the resulting shapefile is empty.

What am I missing?

(I must miss something, because it seems impossible to me that such a basic operation does not work in such a great software, and furthermore this feature has already been mentioned in the 1.8. doku: "Variations in buffering: buffer distance or buffer size can vary according to numerical values provided in the vector layer attribute table for each feature." http://docs.qgis.org/1.8/de/docs/gentle_gis_introduction/09_vector_spatial_analysis_buffers.html)

My project CRS is EPSG:25832, point shapefile the same. QGIS version is 2.8.7LTR (but same symptom in 2.14.3.)

This is what my points look like:

This is how I configure the buffer tool (German GUI):

And this is the result (None):


texts: What open-source web mapping solution I need to use for this situation?

texts:
I will try to explain what I need to do. I have about thirty KML polygons. I need to put them on some satellite map. On the top of these layers (KMLs, satellite) there must go layer with roads. I don't have any "roads" source. I don't care in which format are roads, my primary goal that these roads are on the top of all layers.

My question is what is the best opensource web mapping solution for this kind of a project?
What "roads" layer would you use? Is it possible maybe to put Google ROADS map on the top of satellite imagery and KMLs? If it's not, can I use some vector and where I can find it? And how to put it on the top of other layers?

Please give me also some coding instructions.

Thank you for any help.

EDIT

This is the coding I wrote combining OpenLayers and Google Maps styled maps:

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">

    <meta name="viewport" content="width=device-width; initial-scale=1.0; maximum-scale=1.0; user-scalable=0;" />
    <meta name="apple-mobile-web-app-capable" content="yes" />
    <title>Google styled map &amp; Openlayers</title>

Google styled map &amp; Openlayers

I can see wms layer and KML, but I can't see roads layer. What am I doing wrong?


texts:
I have a few hundred points with ellipsoidal height taken a few years ago. Based on these ellipsoidal heights I had deduced orthometric heights for the same.

Is there a way to find the ellipsoidal height of another point relative to these hundred points on that same ellipsoid? Maybe reconstruct the ellipsoid geometry keeping these hundred as fixed.


texts: UK: Looking to plot house/building numbers on a map, but can't find a data source

texts:
I'm trying to plot UK house numbers on a map but I'm having trouble finding a good data source.

I have download the Ordnance Survey Master Map and the CartographicText layer does contain some house numbers, but nowhere near all of them. I also did a histogram analysis of the data provided, and given the number "12" is more popular than the number "3" and "4" I don't think this can be trusted.

I've also looked into using the Google reverse geocoding API as I have accurate lat/longs for each building, but the results provided are not accurate enough. For example I live at house number 7, but reverse geocoding my position when at home returns house number 5.

Is there any dataset or API which provides accurate house and building numbers for the UK based on a lat/long for a building?


texts: Making street names visible on OSMM WMS, increasing resolution of WMS map at constant scale (Scale Dependent OSMM -SV)

texts: I have my map at a scale where the OSMM background map is not showing the street names. Is there a way to make the street names visible without changing the scale of my map (can't increase scale as showing a city boundary, and cannot increase size of page as needs to be printed off at A2)?


texts:
I am pretty new with all GIS things and I have to prepare homework in QGIS.

I have a shapefile which contains many buildings with their own grades. I want to make an automatic tool, solution, etc. to reach the related pdf reports of each building by clicking.
I read something about actions but I am so confused.

Is there clear way to do this or any source to find about this topic?


texts: get geotiff info like width, height, pixel size with Geoserver REST API

texts: I want read information for certain layer (not tile) like the ones read from gdalinfo. How can I get these info with the rest api?


texts: Emerging hotspot analysis (space-time patterns) in QGIS

texts:
I worked with crime data a while ago and was really impressed by the possibility in ArcGIS to detect significant hotspots and their development over time (Emerging hotspot analysis).

Does anybody know if there is a similar function in QGIS or a plugin?


texts: Trying to read string-based binary (?) code into something like WKT to be read by a spatial database

texts:
I've recently received some parcel data that has "shape" columns which are strings by default and look something like this:

"0F00000060E501A0DFE841410080FF7FD7982741206802C02DE94141808513E02A9A274101000000060000000000000060E501A0DFE841410063FEFFDF982741A031FD7F03E94141808513E02A9A2741206802C02DE94141004C0080199A2741A00B03E028E941410080FF7FD798274140200180E5E841418041FEFFE098274160E501A0DFE841410063FEFFDF98274100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000FFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFEFFF"

What kind of code is this and is there any way to convert it back into something readable by qgis or postgis?


texts: How to create jpg image of map in Geomedia using vb.net?

texts: I am using Geomedia. I want create jpg image of the map. How can it be done using vb.net?


texts: Raster distortion from parallel kriging in R

texts:
I am relatively new to parallel computing in R. I am using a high performance cluster to perform kriging for soil data points. I found this previous post which I used as a template:

How to achieve parallel Kriging in R to speed up the process?

This is the code I have been using, non-parallel, it works fine a produces kriged rasters. For clarity, I am using Random Forests to initially create the spatial predictions from the point data, then run kirging on the residuals.

qs_s5: spatial points data frame of the point data

p2: spatial pixels data frame of the covariate rasters

p <- as(p, Class= "data.frame")
p1 <- p[complete.cases(p),]
p2 <- na.omit(p1)
colnames(p2)[colnames(p2) == "s1"] <- "lat"
colnames(p2)[colnames(p2) == "s2"] <- "long"
v.rf <- variogram(qs_s5$resid.rf~1, data = qs_s5, width = 1000,cutoff=10000)
vm <- vgm(psill = 10, model = "Sph", range =4000, nugget = 20)
vmf.rf <- fit.variogram(v.rf, model = vm)
rf.s_ppm <- krige(resid.rf ~ 1, locations = qs_s5, newdata = p2, model = vmf.rf, debug.level = -1)
gridded(p2)<- ~lat+long
proj4string(p2) <- CRS(utm14)
rf.s_ppm$trend <- p.rf

rf.s_ppm$trend <- p.rf

rf.s_ppm$pred <- rf.s_ppm$trend + rf.s_ppm$var1.pred

rf.s_ppm$lower <- rf.s_ppm$pred - 1.645 * sqrt(rf.s_ppm$var1.var)
rf.s_ppm$upper <- rf.s_ppm$pred + 1.645 * sqrt(rf.s_ppm$var1.var)

writeGDAL(rf.s_ppm["pred"], fname = "./out10/s_ppm_parallel.tif", drivername = "GTiff", type = "Float32")

Here is the code with parallel integrated:

v.rf <- variogram(qs_s5$resid.rf~1, data = qs_s5, width = 1000,cutoff=10000)
vm <- vgm(psill = 10, model = "Sph", range =4000, nugget = 20)
vmf.rf <- fit.variogram(v.rf, model = vm)

gridded(p2)<- ~lat+long
proj4string(p2) <- CRS(utm14)

numCores <- detectCores() -1
numCores

cl <- makeCluster(numCores)

parts <- split(x = 1:length(p2), f = 1:numCores)

clusterExport(cl = cl, varlist = c("qs_s5", "p2", "vmf.rf", "parts"), envir = .GlobalEnv)
clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))

parallelX <- parLapply(cl = cl, X = 1:numCores, fun = function(x) krige(resid.rf ~ 1, locations = qs_s5, newdata = p2[parts[[x]],], model = vmf.rf, debug.level = -1))

stopCluster(cl)

mergeParallelX <- maptools::spRbind(parallelX[[1]], parallelX[[2]])
for (j in 3:length(parallelX)) {
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[j]])
}

mergeParallelX <- SpatialPixelsDataFrame(points = mergeParallelX, data = mergeParallelX@data)

rf.s_ppm <- mergeParallelX

rf.s_ppm$trend <- p.rf

rf.s_ppm$pred <- rf.s_ppm$trend + rf.s_ppm$var1.pred

rf.s_ppm$lower <- rf.s_ppm$pred - 1.645 * sqrt(rf.s_ppm$var1.var)
rf.s_ppm$upper <- rf.s_ppm$pred + 1.645 * sqrt(rf.s_ppm$var1.var)

writeGDAL(rf.s_ppm["pred"], fname = "./out10/s_ppm_parallel.tif", drivername = "GTiff", type = "Float32")

Sdata<-raster("./out10/s_ppm_parallel.tif")
plot(Sdata)

When the raster is generated, the image is distorted to swirls. I understand that this must be do to the way the analysis is cut, distributed, and merged; but I dont know what in the code is the cause nor the solution.


texts: Editing PostGIS database using ArcGIS for Desktop Add-In/Extension without ArcSDE?

texts:
I need to improve an existing solution to be able to do multiuser editing within a PostGIS database from ArcMap.
Available for the solution is ArcGIS for Desktop 10.3 [Basic]. (No non-free Extensions or ArcGIS for Server)

There is already an existing solution implemented as an ArcMap-AddIn. It uses an editable exchange Format (Shapefile) which is read out to the database and updated from the database if there is a change in the Shapefile or PostGIS database.
This solution is prone to errors and needs to be improved.

I know it is possible to read from PostGIS databases in the new ArcGIS for Desktop versions using a query layer. But then the layer can't be edited.
The solution needs to be implemented as an Add-In or an Extension. Whatever is required.

As the “ST-Links SpatialKit” is capable of such a connection. That's why I thought it might be possible to write a ArcMap Extension. Any Advice?

I know similar questions have been asked. But the questions and the answers I found, did relate to old versions of ArcGIS. There has been a lot of change relating PostGIS databases in the resent ArcGIS versions.

I found a clue to a good solution in this answer to a different question.
It links to a coustom implementation of a "Plug-in Data Source" (GDAL-Source). It is a good example I could use as a template for my PostGIS implementation. But this solution didn't implement write support.
The same answer claims, write support for "Plug-in Data Sources" was introdused with version 10.1 SP1. My search didn't find any documentation for this. I can only find the old "read-only" ESRI guidelines for "Plug-in Data Sources".
In the lower part of the comments to the answer, Ragi Yaser Burhum and Petr Krebs try to implement write support. Sadly they didn't resolve there issues (at least publicly). I think this is a branch of their work.
Maybe anybody of you could point me to documentation (or a solution) for write support in a "Plug-in Data Source".