texts: Adding Image Service (i.e., not Map Server) to QGIS 2.18.12

texts:
Novice warning...
I have searched the answers for adding Map Servers to QGIS, but it appears that the link to my service ends with "ImageServer," rather than "MapServer." I cannot convince QGIS 2.18.12 to load any layers from the URL. This is the URL:

http://services.kansasgis.org/arcgis7/rest/services/IMAGERY_STATEWIDE/FSA_NAIP_2014_Color/ImageServer

When I attempt to add this through Layer>Add Layer>Add ArcGIS MapServer Layer, I get this lovely response:

Is there another Add Layer method that would work for this? I have tried all of the other web-based add layer formats with no luck.


texts: Display TMS layer in an ArcGIS Online Map

texts:
When I was using Leaflet to display a tiled map images your answer was very useful,

This was what you answer that time: I only had to add the (-) and everything works good.

Just provide the "url" option instead of tileUrlFunction like:
url: 'https://quadridcmmaps.blob.core.windows.net/tynset/{z}/{x}/{-y}.png'

I am using ArcGIS Online to add a Web Tiled Layer, and I am giving this address, but the images aren't displayed correctly, just like before, please help:

http://localhost/Mapas/2014/bing_roads_en/{level}/{col}/{row}.png

This is the result:


texts: Is it possible to write from a vector memory layer to a Spatialite database?

texts:
Using PyQGIS and QGIS version 2.18.11

I have the following code:

def updateDestination(self, importLayer):
    print "updateDestination called"
    masterName = importLayer.masterLayer
    masterLayer = Utils.getLayerByName(masterName)

    (filePath, fileName) = os.path.split(masterLayer.dataProvider().dataSourceUri())
    filePath = Utils.processTempFilePath(filePath)
    print "Merging into path", filePath
    print "fileName is", fileName

    fileName = fileName[:fileName.index("'")]
    print "fileName is", fileName

    temp = Utils.getLayerByName('Merged')
    print temp

    output = filePath + "/" + fileName
    print "Writing to ", output

    crs = self.iface.mapCanvas().mapRenderer().destinationCrs()
    print crs

    QgsVectorFileWriter.writeAsVectorFormat(temp, output, 'utf-8', crs, 'SQLite', False, None, ['SPATIALITE=YES'])
    QgsMapLayerRegistry.instance().removeMapLayers([temp.id()] )

I already have an existing Spatialite file containing a table with the same name as the file. Let's call it file1.sqlite.

If I change the output to save to "file2.sqlite" it writes perfectly, giving me a new spatialite. I can then import that as a new layer.

My objective is to have it write to file1.sqlite, and override the existing table inside that file with one of the same name, containing the new data (It's the same columns, but has had new rows appended onto the bottom).

When I attempt this, I get no error and no complaints, but it doesn't write. Is there some way to force this, or perhaps a viable workaround?


texts: Generating circular histogram plot of joint directions from vector layer in QGIS?

texts:
I would like to extract the joint direction from a vector layer. Joints are best described 'relaxation fractures' found on rock outcrops are are easily visible if no vegetation exists. The orientation of joints can show the local stress directions caused by faulting or other tectonic activity.

I have identified an outcrop with visible joints and marked the dominant joints as a line within a vector layer - see image. Read lines are the joint directions.. I can see the information I need in the 'Identify Features' option but cannot extract the X1,y1 and x2,y2 information to calculate azimuths.

Are there any plugins to extract azimuths of lines in a vector shapefile within a specified polygon (I would like to extract information within specified areas rather than creating many different 'joint' shapefiles).


texts: Creating a text report detailing attributes on several GIS layers within a polygon

texts: I am trying to create a text report listed by area (taken from a polygon layer) that details any features from other layers (points, lines and polygons) that fall within the listed area.  I have QGIS and FME but I cannot find anything that creates what I need.  Do any of you have any experience of creating this type of output please?


texts:
Unfortunately I have clicked on the edit button in the attributes table and the spreadsheet with the polygons and their attributes has disappeared (see picture).
It is likely that I am making a mountain out of a molehill but I really do not know how to make the spreadsheet reappear.

Thank you in advance


texts:
Is there a way to color a postgis raster using an SLD file?

For example, I want a node server to pull an SLD file from a geoserver instance. Then I want the node server to run a postgis query to retrieve a raster and save it as a png which I've done, but it's gray scale. Is there a way to feed the sld into the query somehow to get color?

I could possibly write a parser that converts the sld into the format needed by ST_ColorMap, but I'm hoping there is an easier way.


texts:
So I have been trying to figure this one out for ages but I just can't seem to do it.

I have 2 rasters. One representing terrain and the other representing an interpolated surface which is beneath the ground.

The interpolated layer is greater than ground level in some places, which in reality is impossible. I need to make the interpolated layer equal to the ground surface layer where it has a greater elevation than the ground surface.

It looks like this in profile view

How can I do this..?


texts: gdal2tiles.py using param: -nb_processes does not use multicore, only one core fully working

texts:
I'm using gdal2tiles to tile huge geoTIFF image, but it takes a long time. I tried using the nb_processes param but only one core is working hard.

My command line:

gdal2tiles.py -e -p geodetic -nb_processes=4 -k -v -w google image.tif TILES/

How can I achieve multicore tiling?


texts: Display GeoJSON on Leaflet web map with GeoDjango

texts:
I have to display my geometry data in my Leaflet based web map and I use GeoDjango to display it, I have loaded GeoJSON data and front map but I can't put my geometry data on my web map.

How can we solve it?

views.py:

from django.shortcuts import render
from django.views.generic import TemplateView
from django.core.serializers import serialize
from django.http import HttpResponse
from .models import jobs

class hompage(TemplateView):
    template_name='index.html'

def job_details(request):
    job=serialize('geojson',jobs.objects.all())
    return HttpResponse(job,content_type='json')

urls.py:

from django.conf.urls import include,url
from .views import hompage,job_details

urlpatterns=[
    url(r'^$',hompage.as_view(),name='home'),
    url(r'^work/$',job_details,name='works'),
    ]

Template.html:

<!DOCTYPE html>

{% load static %}
{% load leaflet_tags %}

{%leaflet_js%}
{%leaflet_css%}
    <title>HARISOHMNAATH</title>

harisohmnaath

{%leaflet_map "gis" callback="window.our_layers" %}


texts: Linking "value maps" together in PostGIS/QGIS where the values are user editable?

texts:
This is related to Limiting available values in dropdown field of QGIS form according to value selected in previous field?. Even though this has been marked as resolved the answer is not what I am looking for as this point towards creating a relation in the project setting.

But, I want to create a relation in 1, or more, database tables to a shapefile itself.

How can I link two, or three, value maps so that the drop down selection of the next field is limited by the value in the previous field? The choices are in a list bound to a database.

The reason why I need it in a choice table list instead of the shapefile itself is because the choice table needs to be updated.

Explanation

For example I have a list to import from excel that has some value that are repeated, but they are distinct because of their subtypes, as illustrated below.

Sometimes there is a subtype, and sometimes not (assumed is 'NULL').

Example:

Begeleidingssnoei - achterstallige snoei
Begeleidingssnoei - bemantelingssnoei
Begeleidingssnoei - beschadigde takken
Onderhoudssnoei - innemen
Onderhoudssnoei - probleemtakken
Onderhoudssnoei - regressie

Because the list is for inventory fieldwork it will not be practical to have these values combined in a single list.

So I want to split it so the value map setting sets the drop down choice list in the next field.

But, importantly these value/s in this next field are different for each previous value.


texts: LineString VS Multi-LineString Performance PostGIS

texts:
Our road dataset (14mil segments, 300,000+ km length) is in Multi-LineString, I am thinking about shifting it all to LineString since all the roads in our data are single geometries and multi part geometries are errors. LineString also provides several different functions that I am really fond of using like ST_StartPoint and ST_EndPoint.

Is there any difference in the performance of LineString and Multi-LineString? and what is the recommended format for storing a large road dataset in PostGIS?

In response to the Comments:
It's not a single Multi-LineString, each one is a different feature. Number of geometries in each feature is 1.
I am not looking for the performance of a specific query.

What I want to know is the difference in speed of full table scans for LineString, and for the same geometry, stored as MultiLineString. Does the geometry type have any impact? The results of explain and the queries I tried, kept on fluctuating and there wasn't any noticeable difference


texts: How to create a polygon based on view extent in QGIS 1.8 ?

texts:
I am looking for a very simple thing. I want to create a polygon based on the current view extent of my project.

In fact I want to

1) open google satellite layers
2) zoom to my area of interest
3) create a polygon based on the 'view extent'

of course, the GE layer is not important, it should work with any view.


texts: Is the ObjectId value guaranteed to be unique after features are deleted and re-created in an enterprise geodatabase?

texts:
I'm working with an enterprise MS SQL Server geodatabase which is accessed via SDE. ArcGIS manages ObjectId values such that they are uniquely assigned as features are added to the geodatabase.
My question is around what happens if a feature (say with ObjectId: 1234) is deleted. Is it a certainty that there will never be a new feature generated with ObjectID: 1234? ie, are previously-assigned ObjectId values ever re-used?


texts: How to make two side by side maps match up in layout view in MapInfo?

texts:
I am working with two maps that have different thematic elements on each- in this case, I made two different columns to display the concentration of chemical found on that location. I display the current standards (2-20, 20-70 and >70) in X, Y, and Z colors in one map, then the other map has the proposed standards which are stricter (environmental regulation goals) as 2-4, 4-10, and >10 as the same colors in the first range. The lowest range is X color, middle range is Y color, and highest range is Z color. I want to make the map images match up exactly, so that people can compare side by side, same scale, same everything, so I can present the differences between the two.

How do I make the two map images match exactly?

I am very inexperienced in MapInfo, my only experience has been mapping and cleaning up data in the tables. I only just started playing with the layout tab about 4 business days ago. I will need a bit more guidance than "just try this", I will have to do a lot of additional research and googling and trial and error to find a solution.


texts: How to calculate a standard deviation across multiple raster layers?

texts:
I have 15 raster layers. Each one shows accumulated precipitation (rainfalls) for the month of April for the same US state during a given year. So, I cover the period from 1997 to 2012. I would like to have an idea of the variability of precipitation during the last 15 years.

I would like for example to calculate the standard deviation at each pixel across the 15 layers. I am guessing I can do it through the raster calculator using sums, subtractions and divisions. But the final formula would pretty long.
Is there a more efficient way to achieve my goal? (So far, I am only a QGis user and I don't know Python...but I can start learning if necessary).


texts: Superimposing high resolution but partial coverage images on a smaller resolution images for use in GeoServer

texts:
For the same area, I have 2 sets of images, the first is high resolution but it doesn't cover fully the area and the images have somehow uneven borders just like this picture.

The second one has smaller resolution, but practically covers my entire area of interest.

I want to superimpose the high resolution images on the smaller resolution image set, but without having these big black images appearing on the map, so I will have to make the black transparent.
I can batch replace with ImageMagick the black color with transparent, but the problem is cannot make PNG converted from the original JP2 to work on GeoServer, while somehow I had managed to make the GeoServer work for the ECW format converted from the original JP2 files. The original .jp2 files were accompanied by .j2w georeference files.
I have to mention that in order to convert the original JP2 files to PNG I use GDAL with a command line gdal_translate -of PNG "C:\Users\userDesktop\file.jp2" "C:\Users\userDesktop\file.png"

I have tried using the png.aux.xml and tried not to, tried to copy from the .prj files from another set that worked and tried not to, all to no avail...

Edit:
Finally with the help of gdalinfo to make sure during each step the georeference was not lost and after various tries, the correct steps are the following
1) Jp2 to png conversion using GDAL
2) png to png_transparent conversion using ImageMagick
3) png_transparent to ecw conversion using GDAL again.


texts:
I want to use ArcPy (for ArcGIS 10.7) and want to change a polyline's point to another coordinates.
I know how to change it (move the point), but when I write like 56.3654275234542375 to set the X value, in ArcPy, the print() function only shows this X coordinate to 9 decimal places.
My questions are:

Can the shapefile store the coordinates longer than 9 decimal places?
How can I change the display length in ArcPy (or ArcGIS)?


texts: Trying to display Sentinel-2 imagery in Google Earth Engine, code working but displaying blank figure

texts:
I'm trying to display SENTINEL-2 imagery in Google Earth Engine, code working but displaying blank figure:

var PuntoQuellaveco1 = /* color: #d63000 */ee.Geometry.Point([-70.477294921875, -17.287709050621917]);
Map.addLayer(PuntoQuellaveco1);

var image = ee.ImageCollection('COPERNICUS/S2')
.filterBounds(PuntoQuellaveco1)
.filterDate('2017-06-01','2017-08-31')
.sort('CLOUD_COVER', false);

var conteo = image.size();
print('conteo de imágenes', conteo);

var mejorImagen = ee.Image(image.sort('CLOUD_COVER').first());

print('La primera con menos nubosidad', mejorImagen);

Map.addLayer(mejorImagen, {bands: ['B4', 'B3', 'B2'], max: 0.3}, 'image');

var fechaAdquirida = mejorImagen.get('DATE_ACQUIRED');

print('Fecha adquirida', fechaAdquirida);


texts: Dissolve polygons in QGIS (GEOS) with tolerance

texts:
I have a polygon vector layer representing the area inundated in a flood. The layer has been produced in a such a way that some of the areas are composed of multiple smaller touching polygons (see screenshot below, labelled "original data").

I would like to dissolve the polygons where they touch. If I use the dissolve geoprocessing tool in QGIS's fTools, which uses the GEOS library, I get some unexpected results (see "dissolve with QGIS fTools"). It seems that GEOS doesn't use a threshold in determining if polygons are touching, as mentioned in these issues:

http://trac.osgeo.org/qgis/ticket/3126

http://trac.osgeo.org/qgis/ticket/3184

http://hub.qgis.org/issues/2806

The same approach in ArcGIS gives the expected result (see "dissolve with ArcGIS 10.1").

Apparently this is a "feature" of GEOS, rather than a bug. Is there any way to fix the data I have to produce the expected result? A minor loss of precision (a few cm?) could be acceptable for me. Ideally the solution would be something I could write into a QGIS plugin, as a pre-processing step.

Related question (with no answers): How to set tolerance value of Geoprocessing tool in QGIS?

Link to sample data as a shapefile: https://www.dropbox.com/s/c8ns7wwfv289ad0/sample.zip


texts:
Anyone have a Java based example of how to implement an asynchronous Geoserver web processing service (WPS)? Based on some presentations that I found, async GeoServer WPS support was added in 2.1.3. I've searched the internet and also the GeoServer and GeoTools websites with no luck.

My environment = PostgreSQL 9.2.7 + PostGIS + GeoServer 2.4.5 + JDK 1.7


texts: Removing smaller area which is less than 10000 inside vector?

texts:
I want to remove smaller area which is less than 10000 inside the vector.But the selected areas(smaller area) only remain after using Data Management>>Eliminate tool in ArcGIS10.1.

Which tool should be used?


texts: Create Polyline feature in Z M Enabled Featureclass

texts:
I want to create a polyline feature in Z M enabled featureclass.

Point1.x= 1000
Point1.Y = 1100
Point1.Z = 1200

PointCollection.addpoint Point1
Point2.x= 5000
Point2.Y = 5100
Point2.Z = 5200
PointCollection.addpoint Point2

But when I try to store feature as Feature.shape = Pointcollection it is giving this error: Item not found in the collection.


texts: What is the difference between qgis-ltr-bin.exe and qgis-ltr-bin-g7.exe for QGIS 3.16?

texts: QGIS 3.16 has two executable files: 'qgis-ltr-bin.exe' and 'qgis-ltr-bin-g7.exe'. Which of these files is the preferred one for running QGIS? Is there any difference between these two files?


texts: Using a variable for map.setcenter in mapbox gl js

texts:
I'm trying to set the center of my map using coordinates passed into the page via the url. I can retrieve the coords, but but get this error when I try to use them:

Uncaught Error: LngLatLike argument must be specified as a LngLat instance, an object {lng: , lat: }, or an array of [, ]

I get the coords like this

function getURLParameter(name) {
return decodeURIComponent((new RegExp('[?|&]' + name +
'=' + '([^&;]+?)(&|#|;|$)').exec(location.search) || [null, ''])[1].replace(/\+/g, '%20')) || null;
}
var center =getURLParameter('location');

open my map with

    mapboxgl.accessToken = 'pk.eyJ1IjoidHJveWIiLCJhIjoiY2pkb3F0cmJnMHA0YTJ4cXB3Mjc0Y3g0eCJ9.3Myq7AzC4yxR8hi9Ubchdw';
    var map = new mapboxgl.Map({
        container: 'map',
        style: 'mapbox://styles/mapbox/streets-v9',
        center: [-88.75, 42.15],
        zoom: 8
    });

and load my layers and set the center of my map in

map.on('load', function() {

map.setCenter(center);

my url looks like this:
http://mywebsite.com/Map.html?location=-84.83762741088867,43.179770125276974&name=otherdata

I can get the name value and use it to specify my map layer.


texts: Assign country names to long and lat in R

texts:
I have gridded data on temperature, from which I want to extract data for specific countries. Put it in another way, I want to kind of assign country names to the coordinates (long and lat). The way I am trying to accomplish it:

library(gapminder)
df = map_data("world") %>% select(long, lat, region)
r = rasterFromXYZ(df)

If I got a raster from df, then I could adjust the resolution and extent to merge it with my data on temperature. Then, I could convert it to a dataframe and subset by country. But when I try to raster df, I get an error message:

"Error in rasterFromXYZ(df) : x cell sizes are not regular".

How can I make it regular? Or is there any other easy way to assign country names to the coordinates?


texts:
How can we add category based styling to geometry point
my application show geoserver rendered map, and all the map POI are directly fetched from database and shown to the map.i need know how we can add different style for different points
please help me

database used is POSTGIS


texts: How to sample points from temporal layerstacking image including null or masked value in Google Earth Engine?

texts:
I'd like to sample points from temporal Sentinel-2 images covering a state, about 13 month's mosaic image were stacked into a new image. When I sampled points on the image, I got nothing. I think maybe some null or masked points led this case. I wonder is there a good solution?

  //Load state borders from assets

var studyArea = ee.Feature(ee.FeatureCollection('ft:17aT9Ud-YnGiXdXEJUyycH2ocUqreOeKGbzCkUw')
    .filter(ee.Filter.eq('id', 'MO'))
    .first()).getInfo().geometry;
var start = ee.Date('2017-01-01');
var finish = ee.Date('2017-12-31');

  // Map the function over one year of data and take the median. // Load
  Sentinel-2 TOA reflectance data.

var s2col = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate(start, finish)
                  .filterBounds(studyArea)
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30));

  // Difference in days between start and finish

var diff = finish.difference(start, 'day')

print(diff)

var temporalResolution = 30 // days

  // Make a list of all dates

var range = ee.List.sequence(0, diff.subtract(1), temporalResolution).map(function(day){return start.advance(day,'day')})

print('range',range)

  // Function for iteraton over the range of dates

var day_mosaics_s2 = function(date, newlist)
{
   date = ee.Date(date)
  newlist = ee.List(newlist)
  var filtered = s2col.filterDate(date, date.advance(temporalResolution, 'day'))
  var image = ee.Image(filtered.mosaic().clip(studyArea)
 .select('B2','B3','B4','B5','B6','B7','B8','B8A','B11','B12'));
  return ee.List(ee.Algorithms.If(filtered.size(), newlist.add(image), newlist))
}

var news2col = ee.ImageCollection(ee.List(range.iterate(day_mosaics_s2, ee.List([]))));

print('mosaicImageCollection_s2',news2col);

var stackCollection = function(collection)
{
  var first = ee.Image(collection.first()).select([]);
  var appendBands = function(image, previous)
 {
    return ee.Image(previous).addBands(image);
  };
  return ee.Image(collection.iterate(appendBands, first));
};

var stacked_s2 = stackCollection(news2col);

print('stacked_composite bands', stacked_s2.bandNames());

var training = stacked_s2.sample(
  {
  numPixels: 500,
  scale:10,
  seed: 0
  }
  );

  print('samples',training)


texts: Geoserver SLD style PointSymbolizer ExternalGraphic issue

texts:
Any idea why this isn't working? It's taken straight from here: http://geoserver.readthedocs.io/en/latest/_downloads/point_pointasgraphic.sld refernced from this page http://geoserver.readthedocs.io/en/latest/styling/sld-cookbook/points.html

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
    xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
    xmlns="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>Point as graphic</Name>
    <UserStyle>
      <Title>GeoServer SLD Cook Book: Point as graphic</Title>
      <FeatureTypeStyle>
        <Rule>
          <PointSymbolizer>
            <Graphic>
              <ExternalGraphic>
                <OnlineResource
                  xlink:type="simple"
                  xlink:href="see-saw.png" />
                <Format>image/png</Format>
              </ExternalGraphic>
              <Size>32</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

My graphic is located in the 'styles' directory of the geoserver data directory, but instead of the graphic I get a grey square. Geoserver version 2.7.1.1


texts: Splitting text from a string field before each delimiter to new fields

texts:
I'm trying to extract text data from a string field containing text delimited by dots and placing it in new fields using the field calculator.

Im using this python function (Taken from How to extract text before a / in QGIS?):

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def func(value1, feature, parent):
    return value1.split('.')[0]

This returns all text before the first dot. Now I'm wondering how to write the function for it to place each delimited text row in a separate field.

Before:

After:


texts: select by location point in polygon with sql query

texts:
I have 2 tables

table a this one a polygon layer with 10000 records

table b this one a point layer with 600000 records

in postgresql geodatabase i select a polygon feature with the query by input x and y

then i selected the one feature polygon i need to select point layer records that

intersect by the polygon feature but this query time 40 seconds and very slow please help me .

select shape from sde.table_b where  st_contains((select shape from st_point(235116.084 ,
2225686.917,10100) as tt , table_a   where st_contains (table_a.shape,tt)='t'),shape)='t'

point data is sample not really very


texts: Densifying many shapes at same time using ArcPy?

texts:
I have all these shapes to densify

I thought it will be better if I code with ArcPy in order to densify all the shapes once at a time.

So I thought of this code:

#creat the list of shapes to densify
liste_ligne =["E1_H20A_H35_D21E", "E1_I_H16A_H18B","E1_K_D16A_D18B","E1_L_H14_H16B","E1_H2_D70_D40C","G1_D50_D55_H40C","E1_H1_D60_D65_H70","E1_M_D12_H12","E1_Bbis_H45","E1_C_H18A","E1_C_H40","E1_D20A_D21A_D35","E1_E_D18A","E1_E_H50","E1_Bbis_D20E1","E1_G2_D21C_H60_H65","E1_LM_LoisirMoyen","E1_Bbis_D20E","E1_F2_D40_D45","E1_LL_LoisirLong","E1_ABIS_H20E_H21A","E1_F1_H55_H21C","E1_Lbis_D16B_D14"]
for ligne in liste_ligne:
    arcpy.Densify_edit(ligne, "DISTANCE", "4 Meters")

The code showed me an error so there's something wrong that I can't figure out.


texts: QGIS raster files not aligned due to different extent

texts:
I have a DEM file of building heights calculated from DSM-DTM in QGIS. Based on this DEM file, I generated a mask that indicates the presence of buildings in matlab. However, when I load the building mask back in QGIS, I found that the two files do not appear on the same extent.

Upon checking their properties, I found that while their CRS are the same, their Layer Extent (layer original source projection) are different.

Both of their Layer Spatial Reference System are as follows:
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 +units=m +no_defs

For DEM file, the layer extent is

570000.0000000000000000,180000.0000000000000000 : 585000.0000000000000000,190000.0000000000000000

For the building mask generated in matlab, the layer extent is

Layer Extent (layer original source projection)

0.0000000000000000,-2033.0000000000000000 :
3983.0000000000000000,0.0000000000000000

I think the problem is that when I carried out image processing in matlab, it just reads the matrix values without retaining its spatial information, thus changing the later extent from 570000 back to 0, as can be seen above.

I have tried to approach it from both QGIS and Matlab, however, with no success so far. I am wondering if anyone has solutions to this or experience with ensuring that the spatial extent remains the same.

Both the files and matlab script can be accessed in this link.


texts:
New to using OpenLayers 3
I want to display two layers from GeoServer: An image layer and a vector layer.  So far I have the following:

// Aerial Tile Layer
        var aerialLayer = new ol.layer.Tile({
            title: "Aerial",
            source: new ol.source.TileWMS({
                url: 'http://server01:8080/geoserver/sf/wms',
                params: {
                    'LAYERS': 'sf:TEST01',
                    'STYLES': ''
                },
                serverType: 'geoserver'
            })
        });

        // View
        var view = new ol.View({
            projection: "EPSG:2236",
            units: 'm'
        });

        // Map
        var map = new ol.Map({
            target: 'map',
            renderer: 'canvas',
            layers: [aerialLayer],
            view: view
        });

However, I am getting the following error. (I am using the ol-debug.js to see a better description of the error):

TypeError: projection is undefined.

Somehow it does not like the projection I am using.  I am not sure why it cannot find that projection.


texts: How do I replace cell values of a raster using values of another raster with differing extent?

texts: I have two rasters (raster1 and raster2). The raster1 extent is smaller than raster2. I want to replace the value of raster2 pixels using raster1 pixels. How can I replace the values?


texts: How to automatically check for lidar flight line overlaps?

texts: I have a series of lidar-based raster products created via FUSION.  Since I am working with low density lidar data, many of the raster images show significant striping where lidar returns overlap on the borders of the flight lines (displayed as red lines).  I need to automatically detect and list these striped rasters for quality control purposes.  What methods can I use to automatically check these rasters for the striped effect shown in the attached image?


texts: Summing values in Virtual Layer returns NULL

texts:
In a Virtual layer I'm trying to sum the values (real) on a table ('PESQUISA') but it returns NULL. What am I doing wrong?
This is the code I'm using:
select sum("renda_mensal") as sum
from PESQUISA


texts: Search functionality in Openlayers/Geoserver

texts:
I am using Openlayers,Geoserver & Postgis.... My postgis tables are having columns as id, name, POINT/POLYGON/LINE.  Each table corresponds to a layer that I am showing in layerswitcher in OL. I am getting these layers via WMS from geoserver.

I want to implement search functionality so that whenever a user gives a string(name value of the table) in search box, he should be zoomed to that particular feature on the map.

Can anyone guide me how to achieve this functionality???
Any help is highly appreciated

The idea I got is:

Step 1: User has to check only one layer(the layer in which he want to search) in layerswitcher before searching . But how to get which layer has been selected??

Step 2: After knowing the selected layer , Using php search that corresponding table for the search string in names column , if it is present get the latitude and longitude values. But how to get these values and what if they are polygons/lines???

Step 3: how to send these latitude longitude values to JS so that I can zoom to that point???

I have completed 2 1/2 steps. I am using st_box2d and getting the bounds of the feature as "BOX(minx miny,maxx maxy)"...... where bounds contain BOX(minx miny,maxx maxy)

the code below is the ajax function for fetching bounds from db. Now how to use bounds variable to zoom to that feature.

 <?-- code for displaying map -->
 function searchfunction(form)
        {
            httpObject = getHTTPObject();
            var searchvalue=form.searchvalue.value;
            if (httpObject != null)
            {
                httpObject.open("GET","searchfeature.php?selectedlayer="+selLayer+"&searchvalue="+searchvalue,true);
                httpObject.send(null);
                httpObject.onreadystatechange = function()
                {
                    if(httpObject.readyState == 4)
                    {
                        bounds=httpObject.responseText;
                    }
                }
            }
        }


texts: Invalid data source error in QGIS while reading .gdb file

texts:
ftp://ftpdata.dnr.sc.gov/gisdata/elev/DEM/

Above is the link to the FTP site provided by SCDNR. It is free data provided by the State of South Carolina. I downloaded the Horrydem and attempted to import it into QGIS but I continue to get the invalid data source error message.

All of the files should be the same. If you're able to import it into QGIS, please let me know how you did it. I've attempted both the vector import command and "drag and drop" as I've seen others mention. None have worked.


texts: Using Geoserver WFS POST and format_options=PRJFILEFORMAT?

texts:
I'm using geoserver (2.4.4) to provide users with shapefile downloads. Some of the software the shapefiles are going to be used in requires .prj-files that follows the "ESRI-standard". A bit of searching lead me to a Geoserver ticket (GEOS-4503), and by downloading esri.properties from there and placing it in my %GEOSERVER_DATA_DIR%/user_projections folder I can issue a GET request to:

http://localhost:8080/geoserver/bw/ows?service=WFS&version=1.1.0&request=GetFeature&typeName=ns:layername&outputFormat=SHAPE-ZIP&format_options=PRJFILEFORMAT:ESRI

and i get back a zipfile with (amongst the other files) a .prj-file with the following:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

This is fine, as my data is in EPSG:4326 and the string above is indeed the "ESRI-type" I want.

However, i also need the user to filter their data, so I've figured that I need to use a filter, and thus POST an XML-file like this:

<wfs:GetFeature service="WFS" version="1.0.0"
                xmlns:wfs="http://www.opengis.net/wfs"
                xmlns:ogc="http://www.opengis.net/ogc"
                xmlns:gml="http://www.opengis.net/gml"
                xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xsi:schemaLocation="http://www.opengis.net/wfs
                http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
                <wfs:Query typeName="ns:layername" srsName="EPSG:4326">

            <ogc:Filter>
                <PropertyIsEqualTo>
                    <PropertyName>version</PropertyName>
                    <Literal>56</Literal>
                </PropertyIsEqualTo>
            </ogc:Filter>

                </wfs:Query>
            </wfs:GetFeature>

to an url like:

http://localhost:8080/geoserver/bw/ows?outputFormat=SHAPE-ZIP&format_options=PRJFILEFORMAT:ESRI

I also get a zip-file with the data back, but this time the .prj file has this content:

GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["E
PSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST]
, AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]]

And this is indeed the OGC-prj-format.

(In addition to setting format_options=PRJFILEFORMAT:ESRI I've also tried to check the "Use ESRI WKT format for SHAPE-ZIP generated .prj files" box in the WFS-page in the geoserver GUI).

Why does the override of prj-file-format work on a GET but not on a POST?


texts: Copying projection of existing defined raster file to define projection of feature class in ModelBuilder?

texts:
For arcmap, can i copy the projection of an existing defined raster file to define a projection of feature class?

I'm planning on doing this inside ModelBuilder.


texts: How to add a field to an OSM file in Python?

texts:
I am writing a program that takes in a .osm file and then calculates some new attributes for that feature.

I am using the python libraries gdal and ogr to help deal with the data. So far OGR has proven to be a very useful library for parsing and accessing the OSM file.

I was wondering if it possible to add an additional field to the OSM file using OGR.

There is a .SetField() function, but it only seems to work on fields that already exist. Is there a way I can design my own field and add it as an extra field to the OSM file?

Bonus points if I can keep the original OSM file intact for testing purposes!

Alternatively, it doesn't have to use OGR, any Python library will do.


texts: List the row of dataframe/attributes in combobox

texts:
I am new to QGIS plugin creation. I have created a 3 Combobox like these Combobox, Combobox_1, Combobox_2. I get the data input from Combobox using the lines
layers = QgsProject.instance().layerTreeRoot().children()
self.dlg.comboBox.addItems([layer.name() for layer in layers])

I stuck here itself. Now I need to select the listed layer from Combobox. Then to list the index of the dataframe in Combobox_1, Combobox_2 which from the Combobox.
For example, I get the input as per above image through Combobox then I need to list like these in Combobox_1, Combobox_2:
 Well_no
 Latitude
 Longitude
 Years
 EC_GEN

Any suggestion?


texts: Seeking location of all Universities and Colleges in North America?

texts:
This seems fairly straightforward but I'm having trouble finding the right data. All I need is the name, address, and lat long of all universities and colleges in North America.

I've come across a list of many universities in America here: http://answers.google.com/answers/threadview/id/578065.html but just by name.


texts: Georeferencing images and losing quality in QGIS?

texts:
I have several plans I need to georeference.  They start out like this:

But after georeferencing, the best I can get is this:

Using these settings:

Where did I go wrong?


texts: the min and the max in raster calculator Arcgis 10.1

texts: Can you help me where I can find or how to calculate the minimum and maximum with raster calculator in arcgis 10.1


texts:
I have an Openlayers 3 app and users can retrieve stored routes from a web service and display them on a map. Currently each route is added to the top of the index but I'd like them to be inserted into an existing Layer Group.

I've tried using these to methods so far but to no avail:

var collection = map.getLayerGroup(plannerLayers);
collection.push(loadedRoute);

map.addLayer(loadedRoute); // this is the only method which works

map.getLayerGroup(plannerLayers).addLayer(loadedRoute);

Layer management in general is something I seem to be struggling with and ideally I'd like to just be able to interrogate the map for vector layers and add discrete functionality to them. Is there a way to programmatically add/remove layers to a group or even move between groups?


texts: Dissolving polygons to single part features with Python (shapely and fiona)?

texts:
I used the code from Dissolving polygons based on multiple attributes with Python (shapely and fiona)? to dissolve a shapefile using shapely and fiona. However the output creates multi-part geometries.

Is it possible to convert the output to single part or generate single part features instead?

I'm attempting to create a single chain using open source tools in Linux. I would also like to avoid calling QGIS tools if that's at all possible.

Here's the code sample:

from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
with fiona.open('cb_2013_us_county_20m.shp') as input:
    # preserve the schema of the original shapefile, including the crs
    meta = input.meta
    with fiona.open('dissolve.shp', 'w', **meta) as output:
        # groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'STATEFP' field
        e = sorted(input, key=lambda k: k['properties']['STATEFP'])
        # group by the 'STATEFP' field
        for key, group in itertools.groupby(e, key=lambda x:x['properties']['STATEFP']):
            properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
            # write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
            output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})


texts: When adding addin extension to project Startup function is never hit with ArcObjects and C#?

texts:
I have created an Esri Addin toolbar in VS2010 to add layers and search an address file, this works fine. I want to add an extension that will replace text on a template MXD with the current users details. I have all the logic working in a base class and would like to add an extension that listens for the user to switch to layout view and then calls a function in my base class. I have everything working except that when I right click the project and add an extension the extension code never loads. I am probably missing something pretty basic, could someone please point it out to me.

[Guid("ff8f632e-79c3-4168-a591-441bb8f6c2d7")]
[ClassInterface(ClassInterfaceType.None)]
[ProgId("DCC_Toolbar_10.addUserDetailsextension")]
public class addUserDetailsextension : IExtension
{
    #region COM Registration Function(s)
    [ComRegisterFunction()]
    [ComVisible(false)]
    static void RegisterFunction(Type registerType)
    {
        // Required for ArcGIS Component Category Registrar support
        ArcGISCategoryRegistration(registerType);

        //
        // TODO: Add any COM registration code here
        //

    }

    [ComUnregisterFunction()]
    [ComVisible(false)]
    static void UnregisterFunction(Type registerType)
    {
        // Required for ArcGIS Component Category Registrar support
        ArcGISCategoryUnregistration(registerType);

        //
        // TODO: Add any COM unregistration code here
        //
    }

    #region ArcGIS Component Category Registrar generated code
    /// <summary>
    /// Required method for ArcGIS Component Category registration -
    /// Do not modify the contents of this method with the code editor.
    /// </summary>
    private static void ArcGISCategoryRegistration(Type registerType)
    {
        string regKey = string.Format("HKEY_CLASSES_ROOT\\CLSID\\{{{0}}}", registerType.GUID);
        MxExtension.Register(regKey);
    }
    /// <summary>
    /// Required method for ArcGIS Component Category unregistration -
    /// Do not modify the contents of this method with the code editor.
    /// </summary>
    private static void ArcGISCategoryUnregistration(Type registerType)
    {
        string regKey = string.Format("HKEY_CLASSES_ROOT\\CLSID\\{{{0}}}", registerType.GUID);
        MxExtension.Unregister(regKey);
    }

    #endregion
    #endregion
    private IApplication m_application;

    #region IExtension Members

    /// <summary>
    /// Name of extension. Do not exceed 31 characters
    /// </summary>
    public string Name
    {
        get
        {
            //TODO: Modify string to uniquely identify extension
            return "addUserDetailsextension";
        }
    }

    public void Shutdown()
    {
        //TODO: Clean up resources

        m_application = null;
    }

    public void Startup(ref object initializationData)
    {
        m_application = initializationData as IApplication;
        if (m_application == null)
            return;

        //Listening to document events
        SetUpDocumentEvent(m_application.Document);
    }

    #endregion

     #region Document events
    //Event member variable.
    private ESRI.ArcGIS.ArcMapUI.IDocumentEvents_Event m_docEvents = null;

    //Wiring.
    private void SetUpDocumentEvent(ESRI.ArcGIS.Framework.IDocument myDocument)
    {
        m_docEvents = myDocument as ESRI.ArcGIS.ArcMapUI.IDocumentEvents_Event;
        m_docEvents.NewDocument += new ESRI.ArcGIS.ArcMapUI.IDocumentEvents_NewDocumentEventHandler(OnNewDocument);
        m_docEvents.OpenDocument += new ESRI.ArcGIS.ArcMapUI.IDocumentEvents_OpenDocumentEventHandler(OnOpenDocument);
    }

    //Event handler methods.
    void OnNewDocument()
    {
        System.Diagnostics.Debug.WriteLine("New document event (C# Sample)");
    }

    void OnOpenDocument()
    {
        System.Diagnostics.Debug.WriteLine("Open document event (C# Sample)");
    }
    #endregion
}


texts: Public transport drive times in Europe/ Germany possible with Deutsche Bahn?

texts:
Background:

There are many data sources for calculating drive times for driving, walking or cycling, for example Google Maps DistanceMatrix, ESRIS WebServices (no bicycle at the moment) etc. They all have there advantages and disadvantages but with a little effort you can get the needed information.

But, on the one hand, there seems to be no comprehensive data source for public transport information in Europe and especially in Germany. I have tried some alternatives like Google Maps DistanceMatrix - problem here is that only some parts of Germany support public transport outside the long-distance traffic - or OpenStreetMap - partially incomplete, outdated or wrong information.

On the other hand you have a mighty service from the Deutsche Bahn, like http://www.fahrplanauskunft.de/bin/query.exe/d. This database is not limited to their own "products" like Deutsche Bahn trains but also shows local public transport information or information of their regional competitors. On top of all, their information is not limited to the country border of germany but comprehensively cover large parts of Europe.

Question:

What options are there to use their information in ArcGIS/ ArcPy/ Python/ etc. for calculating drive times from one point to an other?

I found a project which also want to participate on their information (OpenBahn API), but this project has stoped. Are there any comparable APIs or possible opportunities (like json, http get etc.) to use this information without a manual search on one of the websites of Deutsche Bahn? At the beginning my first goal is only to send one information like starting point, ending point and arrival/ departure time and receive a computer readable answer like duration or maybe the frequency of daily connections.


texts:
We need to add esri-tile layer in leaflet .At first ,we tried proj4leaflet with esri-leaflet 2.0.0 and leafelt 1.0.0. It workded fine.

<!DOCTYPE html>

     <link rel="stylesheet" href="https://cdn.jsdelivr.net/leaflet/1.0.0-rc.2/leaflet.css" />

Our program is built on leaflet 0.7.7 only working with esri-leaflet 1.0.4. So we also tried it with leaflet 0.7.7 and esri-leafelt 1.0.4 ,but it complains L.esri.TiledMapLayer is using a non-mercator spatial reference .

Is it possbile to run proj4leaflet without upgrading leafelt and esri-leaflet ？


texts: Seeking plugin for QGIS to create iso-distance polygons longer than 120km

texts:
I'm new on QGIS and I have been using the ORS plugin to create isochrones and iso-distance polygons. However, through Openroute Service I cannot create iso-distances longer than 120km.
Is there any other plugin to overcome the limit of 120km?


texts:
I am attempting to move simple Geoprocessing routines from ESRI based processes to SQL Server. My assumption is that it will be far more efficient. For my initial test I am working on an intersection routine to associate overlapping linear data.

In my WCASING table I have 1610 records. I am trying to associate these Casings with their associated mains. I have ~277,000 Mains. I have ~1,600 Casings.

I am running the query below to get a general sense of how long it will take to find individual matches. This query returned 5 valid intersections in 40 seconds.

SELECT Top 5 [WCASING].[OBJECTID] As CasingOBJECTID,
    [WPUMPPRESSUREMAIN].[OBJECTID] AS MainObjectID, [WCASING].[Shape]
FROM [dbo].[WPUMPPRESSUREMAIN]
    JOIN [WCASING]
        ON [WCASING].[Shape].STIntersects([WPUMPPRESSUREMAIN].[Shape]) = 1

My Primary questions;

Will this process faster depending on the search order?

Finding 'A' inside of 'B' vs
Finding 'B' inside of 'A'
Initial return on 5 records from these datasets is that it does not matter

Will this process faster, if I first buffer to limit to a smaller main set and then search?

Can I use SQL Server Tuning to work with Geometry based queries?

SELECT WCASING.OBJECTID AS CasingOBJECTID,
    WPUMPPRESSUREMAIN.OBJECTID AS MainObjectID, WCASING.UFID AS UFID,
    WPUMPPRESSUREMAIN_IPS.UFID AS MainUFID, WCASING.SHAPE
INTO WCASING_INTDefsV6
FROM WCASING with (index([FDO_ShapeWC]))
    INNER JOIN [WPUMPPRESSUREMAIN_IPS] ON
        [WPUMPPRESSUREMAIN_IPS].Shape.STIntersects(WCASING.SHAPE) = 1

This new query has improved definitions.

Now both tables have spatial indexes
Previously Casing Table (Smaller) did not have a Spatial Index

It did contain a Non-Clustered Index

The query also has the with index statement.

The new query took 37 minutes.
The old query took 44 minutes.

I was hoping for better results and will keep testing.


texts: Carto.js Infowindow Listeners open/close Events

texts:
I have a website with the carto.js library and some layers loaded. If I click on a feature, an Infowindow is opening up.

I know that there is an event for the feature click:

sublayer.on('featureClick', function(e, latlng, pos, data) {
          alert("Hey! You clicked " + data.cartodb_id);
    });
  }).on('error', function() {
        //log the error
});

But I have not found any information on how to listen for the close event. Is it possible to create a listener on the close event of the infowindow?

Usecase: Change map content while the popup is open.


texts: Problems with intersecting two polygons layers in QGIS

texts:
I have two polygons layers, and I'm looking to find the % Area of layer A polygons that are placed inside layer B polygons. (some of them are totally inside ,so I expect to get the result "1", and some of them are half in and half in other polygon, so I expect to get "0.5" for both polygons of layer B)

When i'm using the Intersection tool (Vector->Geoprocessing tools), I get the following error: "Input layer A contains invalid geometries (Feature 1426). Unable to complete intersection algorithm."

When I try to use the "Join Attributes by Location" tool (Vector->Data Management Tool),using "Intersects" Geometric Predicate, I get a table with all layers A that intersect layer B, that's great but I can't calculate the % area of layer A polygons that are inside layer B polygons.

I'm using QGIS 2.18.4. Both layers are CRS EPSG:4326, WGS 84.
Any advice?


texts: Creating drop-down list from dynamic field values in ModelBuilder?

texts:
I am trying to create a drop-down list containing all the values from a particular field, that I can use as an input for the Select by Attributes tool in ModelBuilder.

The answer I’ve seen to similar questions is to hard-code the list using the String Variable.

i.e. Right-click in white space > Create Variable… > String  > Choose to Model Parameter > Model > Model Properties… > make String filter = Value List > enter values one by one > Then use %String% as an in-line variable in a tool expression.

However, as the field is regularly updated by other people, I don’t want to hard-code any values.

Is there any way that I can perhaps, iterate through a layer (called "Boundary") and collect the values from a given field ([Name] field) as a String or Value(?) that can be used as an in-line variable in another tool, with ModelBuilder (or with Python script)?

The outcome being a drop-down list filled with values from the Boundary layer's [Name] field in the Select by Attributes tool.


texts: Is it possible to export a filtered layer in ArcGIS Desktop?

texts:
I'm considering purchasing ArcGIS as an ArcGIS noob. (Their trial is currently limited to Coronavirus support clients).

I would be using the software to generate shapefiles of US zip codes by state. i.e. I will have 50 separate files. I looked at the US Census Bureau's shapefile offerings and they only have a single shapefile available with all US Zip codes. (Actually they are ZCTAs - zip code tabulation areas). The resulting file is very large and slow for our website to handle. That's why we want state by state zip code shapefiles.

I played around with ArcGIS online and see it's possible to import a US Zip codes layer to my map and then filter by state.

Please correct me if I'm mistaken, but I read it is not possible to download shapefiles from ArcGIS online. That being said, I am curious if I will be able to create this same filtered layer map in ArcGIS Desktop and then download the resulting map as a shapefile. If so, please point me toward a tutorial or documentation on how to do so.

Here is an article describing how to do something similar.
And some other potentially relevant information.

At a more detailed level, we are actually first converting our shapefiles to vector tile sets with MapTiler and then importing these tile sets as protocol buffer files into our client. If it's possible to directly convert a filtered layer in ArcGIS into a vector tile set (i.e. skip the shapefile step) that we can then use in our site as a protocol buffer, that would be ideal.


texts: Why are heights given at various points on roads?

texts:
In the OS 1:25000 map below I have circled some numbers next to roads. The key says "Ground survey height" and there is another entry (in brown) for "Air survey height", but since all the data in the map comes from a survey of one kind or another, I'm not really sure why this is relevant. You can already tell the height within 5m from the contours, and you can make a reasonable estimate if you want more precision than 5m.

I have tried to think what they might used be for, perhaps for knowing drainage directions, but it seems too little information for that (i.e. you'd need a more detailed and specialized map for that). So...

What are these specific heights used for (why do they appear on this type of map)?

Is there any reason those particular points are used, they seem fairly arbitrary?

Above map is © Crown copyright (2016) Ordnance Survey. Used for private research in accordance with the principles of fair dealing as set out in the Copyright, Designs and Patents Act 1988.


texts:
I'm trying to import a large GeoTIFF (~1.2 GB) into SAGA GIS for processing as a grid using Import/Export - GDAL/OGR -> GDAL: Import Raster, but it keeps crashing. I'm using the 32bit OSGEO4W SAGA installation, on Windows 7.

This is the only error message I can see before the GUI freezes, it seems to be a problem with memory.

Does anyone know how I can change memory allocation settings in SAGA to avoid this crash? I'm sure the station I'm working on has sufficient memory as I can view this raster fine in QGIS and ArcMap. And I'm sure SAGA is working fine as it imports smaller rasters no trouble.

UPDATE
I have since installed 64 bit SAGA GIS and this imports the raster no problems. I found this resource on the SAGA forums, explains the memory limitations nicely.


texts: Extracting coordinates of polygon vertices in ArcMap?

texts:
I've got about a dozen polygons in a feature class loaded in ArcMap 10, all in geographic WGS 1984.

How do I easily obtain the coordinates associated with each vertex of each polygon in that feature class?

Ideally I would like them to be nicely tabulated in spreadsheet format.


texts: How to symbolize raster with field values in QGIS or GRASS GIS?

texts:
I have a raster with ten field values.

In ArcGIS I can open and select one of the field values for symbology and it show me the values from that field.

How can I handle this in QGIS?

I share a picture for a better comprehension.

ArcGIS

QGIS

Screenshot showing the six files of raster file.


texts:
I'm sub-setting a GeoTIFF with the GDAL Python binding and numpy. Essentially, I am providing a central location (coordinates of a flux tower) from which I produce a subsetted image. I then want to write this to a new GeoTIFF. However, my current code does not correctly georeference. I'm having a hard time digesting the rather sparse and confusing GDAL documentation and so would love for some help. I have read that this might be possible with a gdal_translate or a VRT, but can not figure out any more from the docs.

So, here is how I load the GeoTIFF in question and produce the subset:

datafile = gdal.Open("%s" % filenames_matched[site][filen] )
data = datafile.ReadAsArray()
info = datafile.GetGeoTransform()
# GDAL provides the top-left coordinate (lat, lon)
#Site_lon and site_lat are the coordinates for the fluxtower site.
# Produce a relative coordinate for the fluxtower (eg relative to gdal ul(lat,lon))
lon_rel = site_lon - info[0]
lat_rel = site_lat - info[3]

#Get pixel size (resolution)
res_we = abs(info[1])
res_ns = abs(info[5])

#Now convert the relative position from degrees to pixels..
pixel_lon = float ( lon_rel) / res_we
pixel_lat = float ( lat_rel ) / res_ns

# Will, of course, have to round this to the nearest whole pixel...
pixel_lon = round( pixel_lon )
pixel_lat = round( pixel_lat )

# Now need to outline the actual buffer of 24 on each side...
buffer_lon_east = pixel_lon + 24.5
buffer_lon_west = pixel_lon - 24.5
buffer_lat_north = pixel_lat + 24.5
buffer_lat_south = pixel_lat - 24.5

#Produce subset array
subset = data[buffer_lon_west:buffer_lon_east,buffer_lat_south:buffer_lat_north]

And then save the new subsetted GeoTIFF:

Nx, Ny = subset.shape
driver = gdal.GetDriverByName("GTiff")

projection = datafile.GetProjection()    # Same projection as input (WGS84)
ds = driver.Create("%s" % composited_filename, Nx, Ny, 1, gdal.GDT_Float32)
ds.SetProjection(projection)

# Where I imagine the problem is:
# Sets new coordinates from derived ul lat, lon for the subset
ds.SetGeoTransform([info[0] + lon_rel, info[1], info[2], info[3] + lat_rel, info[4], info[5])

ds.GetRasterBand(1).WriteArray(subset)
ds = None


texts: Making RasterToASCII_conversion observe environment setting for output cell size?

texts:
I have two rasters (an ESRI grid and an IMAGINE image) that differ in projection, extent and cell size. I want to convert the IMAGINE to an ASCII, and want it to have the extent, cell size and projection of the ESRI grid. I've set global environment variables (extent, output CRS, cell size, snap raster) to match the values in the ESRI grid, and the ASCII that results from RasterToASCII_conversion has the correct projection and extent, but ignores the cell size setting.

However, if I use the Raster and ApplyEnvironment tools (see code below) from the Python console, the output is as expected (correct cell size, extent, etc).

import arcpy
arcpy.env.extent = 'c:/path/to/grid'
arcpy.env.cellSize = 'c:/path/to/grid'
arcpy.env.outputCoordinateSystem = 'c:/path/to/grid'
arcpy.env.snapRaster = 'c:/path/to/grid'
r = arcpy.Raster('c:/path/to/imagine')
arcpy.CheckoutExtension('Spatial')
rclip = arcpy.sa.ApplyEnvironment(r)
rclip.save('c:/path/to/outfile')

While I could do the above and then convert the resulting grid to ASCII, I'd prefer to do the whole procedure without using Spatial Analyst. Is this possible?

I'm using ArcGIS 10.


texts: QGIS 3 Dragging layer in table of content with keyboard

texts:
Sometimes it is hard to drag a layer with a mouse.
Is there a way to drag a layer (сhange its position relative to other layers) in a QGIS 3 table of content using only the keyboard?


texts: Getting pixel size without opening raster using Python

texts:
I want to get the pixel size (resolution) of a raster without opening up the raster in gdal. I currently have my script opening the raster in GDAL using gdal.Open() and getting the resolution that way. However this process is slow especially when my TIFFs are large and I have many of them.

from osgeo import gdal
data_path = "some/path/to/raster.tif"
ds = gdal.Open(data_path)
# md = ds.GetMetadata() # read metadata
gt = raster.GetGeoTransform()

pixelSizeX = gt[1]
pixelSizeY =-gt[5]

Is there a way to grab this info from the metadata or something without opening the file?


texts: Plot geometry stored in Postgresql/Postgis with GDAL/OGR

texts:
I have a table in which I store a geometry (point, polygon ect). I want to wrote a python program that "draw" those geometry.

I tried to use pgsql2shp and then gdal_rasterize, but the fact that I have mixed geometry types prevent me to progress any further.
And when I try with only one geometry, gdal_raserize raise this error :

Attempt to create 0x0 datasete is illegal, size must be larger than zeo.

Those are the command I run :

pgsql2shp -f shapefile.shp -h localhost -u jenkins -P jenkins -p 5432 video "SELECT geom FROM state WHERE name LIKE \'A%\'"
gdal_rasterize -b 1 -b 2 -b 3 -burn 255 -burn 0 -burn 0 -l shapefile shapefile.shp work.tif

The main objective, is to draw those geometry on a map.


texts: How should I correctly exaggerate a hillshade raster (ArcGIS10, Spatial Analyst)?

texts:
I once wanted to generalize a high resolution raster and instead of MEAN i accidentally selected the default SUM technique, using a cell factor 5 (Lidar DEM is 1m resolution).  When I used this result as input to the hillshade funciton what I got was what I can only describe a highly exaggerated hillshade.   I deleted the results and recreated the generalized DEM but now I have the need to show more exaggerated terrain characteristics in a topographic map and it seems to produce the results I want albeit generalized.  Is this a common way to do this or should I be using a different technique?
How would I exaggerate a DEM without having to aggregate the raster?

Hillshade created from generalized (Mean, Cell factor 5) LiDAR DEM

Hillshade created from generalized (SUMm Cell factor 5) LiDAR DEM


texts:
How can include a title legend (instead of legend) in geonode?

I have tried with everything but I can't do it.


texts: How can I fix incorrect links in my QGIS install?

texts: Today I installed QGIS 1.8.0 and I have the same error as this one: How can I fix QGIS 1.8 startup issues due to missing qwt5.dll library?. The sollutions how ever don't seem to work. What I have found out about it is that certain files and subfolders are in folders where they are not supposed to be. For example: the QGIS-executable is in c:\program files\quantum gis lisboa\apps\qgis\bin\ while it should be in c:\program files\quantum gis lisboa\bin\ where also the qwt5.dll-file can be found. By copying the files from the first folder to the latter I am able to start QGIS. I still get errors though, for example that plugins can't be found. This again, I can solve be copying folders. Now I can use QGIS more or less, but it doesn't run like it is supposed to. Does anyone know how I can fix this problem?


texts: Using unplaced annotation without editing

texts:
I'm working with a layer that isn't editable, and when attempting to convert labels to annotation, there were some unplaced annotations.
In the past, I would just open the unplaced annotation window and add them. I have discovered that in 10.0, the annotation has been moved into the editing environment.
Since I can't edit the layer, I can't access the unplaced annotation. by the same token, I can't create database stored annotation since the destination is not an editable environment, and can't be changed.
Apart from exporting the layer in question into an editable location, is there a workaround for map-stored annotation that can be accessed outside of an edit session?


texts:
I wanted to install the plugins CLUZ and Marxan in QGIS, but both do not appear when you search for them in the plugins sections...

Does anyone know why this is the case, or how can I use CLUZ then?

(I downloaded the latest version Madeira 3.4.2 today)


texts: TOC layer symbol changed EVENTS in ArcMap

texts:
I want to do something when user changes one of TOC layer's symbol. I am using .NET in arcMap 10.0 . How can it possible with ArcObjects ?

Is there any events that will fired when user change symbol of a layer in TOC?


texts: Importing Data from API-Website into QGIS

texts:
I am trying to use the map data, found on this website

https://fixmyberlin.de/planungen

and import it to QGIS. Specifically I am looking for the bike paths shown there. I have contacted them and they sent me a link to the API version of their website:

https://api.fixmyberlin.de/api/projects?page_size=200

Since I am new to QGIS and that whole topic, I really don't know what to do further.


texts:
So! Really basic question. I want to make a Torque map – easy enough – but I want it to display the various datapoints in different colors depending on the value of annual_irr.

Here's what I thought I'd enter into the cartocss editor to make it work (where annual_irr is a number from -.06 to .24), but this just displays the datapoints with no color at all:

Map {
-torque-frame-count:512;
-torque-animation-duration:30;
-torque-time-attribute:"bldg_bought";
-torque-aggregation-function:"count(cartodb_id)";
-torque-resolution:2;
-torque-data-aggregation:cumulative;
}

#fake_porchlight_thing_for_cartodb_larger_bldgs_1{
  comp-op: lighter;
  marker-fill-opacity: 0.9;
  marker-line-color: #FFF;
  marker-line-width: 1.5;
  marker-line-opacity: 1;
  marker-type: ellipse;
  marker-width: 6;
  [ annual_irr <= 1 ] { marker-fill: #1a9850; }
  [ annual_irr <= 0.1 ] { marker-fill: #fff2cc; }
  [ annual_irr <= 0 ] { marker-fill: #d73027; }
}

Any idea what I'm doing wrong? Why isn't the CartoDB thingy interpreting the annual_irr marker-fill rules correctly?


texts: How to reset Network Dataset in Network Analysis Properties with python (arcpy)?

texts:
I have mxd which is used to upload network service.
Mxd uses our network dataset to do Closest Facility and Routing Analysis. All features of this Network Dataset are listed in the table of content as layers.
This mxd is used to create service definition file for our portal.

I have a script which updates layers in this mxd with the latest road network and uploads a new service definition file. Script works fine.

However I have recently realized that updating the layers in the map does not actually updates the sources for network dataset that is used for the analysis.

Does anyone know how to update source for the Network Dataset with python?
I have inserted below the screenshot of the menu (with button highlighted in red). This is what I am looking to do with arcpy.


texts: Pairing elements of two lists in Google Earth Engine

texts:
In Google Earth Engine, I have two lists:

var list1 = ee.List.sequence(50, 100, 50);
var list2 = ee.List.sequence(2,3,1);

How to make a new list containing every possible pair from list1 to list2 ?

I tried zip() function which gave results as follows:
[
  [
    50,
    2
  ],
  [
    100,
    3
  ]
].

Pairs [50,3] and [100,2] are missing from this new list.
Is there any way out?


texts: How to place items horizontally below heading in Legend in ArcMap?

texts: I'm working with ArcGis 10.2. I want to display a legend which is composed of Title ("Legend"), Heading ("Layer Name") and Items. I want the items (in this case three) to be horizontally disposed under the heading, so I set 3 as number of columns, but this way one column was assigned to the heading and two to the items (resulting in the third item to be in the same column of the heading or another item, as you can see in the picture below) while I want the heading to be above all three columns, each of which assigned to an item. How can I do?


texts:
I am attempting to extract natural break values from a raster using Python.

The intended process goes: Arcpy raster object to NumPy Array (using RasterToNumPyArray) to natural break values (using PySAL Natural Breaks function).

import arcpy, pysal
from pysal.esda.mapclassify import Natural_Breaks as nb
# code to create greenIndex arcpy Raster object here
greenArray = arcpy.RasterToNumPyArray(greenIndex)
breaks = nb(greenArray,k=2,initial=20)

This code returns the error, "ValueError: matrix must be 2-dimensional".

As far as I know, greenArray is a 2-dimensional array.


texts: Convert non-spatial Oracle table to spatial Oracle table in a different database

texts:
I'm trying to use ogr2ogr to convert data from a non-spatial table (has X & Y values) in an Oracle database (1) to a spatial format in a different Oracle database (2).  This will be scheduled to run daily to overwrite the table in Oracle database (2) each time.

I have it working using a ogr2ogr's Virtual File Format, with the file Oracle.vrt which is like:

<OGRVRTDataSource>
  <OGRVRTLayer name="TEST_X_Y">
    <SrcDataSource>OCI:userid_1/password_1@database_instance_1:TEST_X_Y</SrcDataSource>
    <SrcLayer>TEST_X_Y</SrcLayer>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>EPSG:27700</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X_TEXT" y="Y_TEXT"/>
  </OGRVRTLayer>
</OGRVRTDataSource>

Then using the following ogr2ogr command:

ogr2ogr -update -overwrite -f OCI OCI:userid_2/password_2@database_instance_2:TEST_OGR_IN Oracle.vrt TEST_X_Y -lco OVERWRITE=YES -lco SRID=27700 -nln TEST_OGR_IN

Now this works and creates the spatial table in the Oracle database (2) fine, but when I run the ogr2ogr command a second time I get an error like:

ERROR 1: ORA-00001: unique constraint (MDSYS.UNIQUE_LAYERS) violated
ORA-06512: at "MDSYS.SDO_GEOM_TRIG_INS1", line 43
ORA-04088: error during execution of trigger 'MDSYS.SDO_GEOM_TRIG_INS1' in
  INSERT INTO USER_SDO_GEOM_METADATA VALUES ('TEST_X_Y',
...

The table (TEST_OGR_IN) still gets overwritten/updated in the Oracle database (2) but the problem appears to be that ogr2ogr has inserted values into the USER_SDO_GEOM_METADATA table on the Oracle database (1) for the non-spatial table (TEST_X_Y) and is violating an index/constraint?

Although the process works, I'd like to understand why the error is happening and how I can avoid it.


texts: QGIS: What is the purpose of the Browser "Fast Scan this Directory" option?

texts:
At 3.10.5, I notice that there is a Browser checkbox for Fast Scan this Directory.  However, when I check it, nothing seems to change in the way the Browser operates:

The 3.10 documentation does not state its purpose, nor did an internet search turn up anything.

So, what is the function of the Fast Scan option?


texts: Error running MapServer export in Quantum GIS

texts:
I have a project in Quantum GIS that I need to export to a .map file using the plugin MapServer export.  I fill all the forms but I get this error when running the app.  Any ideas of what is happening?

Here is a copy of the error:

Traceback (most recent call last):   File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins\mapserver_export\mapserverexport.py", line 132, in ok_clicked
    self.saveMapFile()   File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins\mapserver_export\mapserverexport.py", line 148, in saveMapFile
    print "Creating exporter using '%s' and '%s'" % (self.dlg.ui.txtQgisFilePath.text(), self.dlg.ui.txtMapFilePath.text()) UnicodeEncodeError: 'ascii' codec can't encode character u'\xed' in position 49: ordinal not in range(128)

Versión de Python:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

Versión de QGIS:
1.8.0-Lisboa Lisboa, 6416f38

Ruta de Python: ['C:/Users/SergioVM/.qgis//python/plugins\\sextante', 'C:\\Users\\SergioVM\\.qgis\\python\\plugins\\elevation', 'C:\\Users\\SergioVM\\.qgis\\python\\plugins\\cswclient', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/SergioVM/.qgis//python', 'C:/Users/SergioVM/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\\Windows\\system32\\python27.zip', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\DLLs', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\qgis\\bin', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', '~/.qgis/python', 'C:\\Users\\SergioVM\\.qgis\\python\\plugins\\openlayersov\\logic', 'C:\\Users\\SergioVM\\.qgis\\python\\plugins\\openlayersov\\gui', '/usr/share/qgis/python', 'C:\\Users\\SergioVM\\.qgis\\python\\plugins\\QuickMultiAttributeEdit/forms', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\qgis\\python\\plugins\\fTools\\tools', 'C:/SIG-SIXAOLA/Cartograf\xc3\xada/base_datos_final/Quantum_GIS_BD']


texts: Why does my PyQGIS code for saving selected features result in altered geometries?

texts:
In a plugin I want to save the selected features in a new shapefile. I search and I try  (How can I create a new layer from the currrent selection? to create a new layer and write the selected fields and features in that:

    vl = QgsVectorLayer("Polygon", "temporary_points", "memory")
    pr = vl.dataProvider()

    # add geometry and fields
    pr.addAttributes(fields)
    pr.addFeatures(features)

But this result is not what I wanted. I get this:

In this image, the selected polygon is yellow and the pink polygon is the saved one as another shapefile. So the geometry is not completely correct and I need the perfect copy of the original.

Anyone knows how to get it?

I can use the algorithm from Processing Toolbox, "Save selected features".


texts:
Does anyone know where I can get or find shapefiles for pipelines (USA location)? Preferably free as I've never done pipelines before and I'm meeting a client next week that is doing pipeline work.

How much different is pipeline from OG well location/placement etc.?  I use ArcEditor 10.0.


texts: GeoServer creates thousands of .tmp files

texts:
I noticed that I currently have 34736 .tmp files in GEOSERVER_ROOT/data directory along with the usual files and directories.
From what it seems, these files are copies of the global.xml file, as their names are:

global.xml.42e17d18-1e2a-4025-91e5-b43cda42d155.tmp
global.xml.cbf68061-bab3-4719-bb9c-bf77040f47ae.tmp
global.xml.25d2c6f8-445e-4854-9fb8-b2902e0b44be.tmp

The only difference between these files and the original global.xml is the value of updateSequence,
which increases by one with every new file created:

Global.xml:
<updateSequence>138667</updateSequence>

.tmp files:
<updateSequence>207176</updateSequence>
<updateSequence>207177</updateSequence>
<updateSequence>207178</updateSequence>

The value from the last created file corresponds to the value listed on Server Status in GeoServer GUI:

Update Sequence 210867

From GeoServer User Manual:

  "Update Sequence  Refers to the number of times (60) the server
  configuration has been modified."

Possible explanation: I have an automated script that fetches some files from a FTP, processes them and uploads them to Geoserver via curl. The script runs every hour.

So the configuration is being modified, but shouldn't the global.xml file be updated?

I could create a .bat file that deletes them once in a while, but I would like to find out why they are created in the first place.

So my question is:

Why does GeoServer create these temp files (and forgets to delete them) and doesn't update the original global.xml?

- EDIT -

I tried to take it step by step to see what happens.
Apparently, if I create a new workspace and I upload a file in the newly created coverageStore, it works without any errors and no global.xml.tmp files.

But if I do that for all the files from a folder..

import urllib
import urllib2
import base64

base64string = base64.encodestring('%s:%s' % ("x", "x"))[:-1]
authorizationHeader =  "Basic %s" % base64string;

baseURL = "http://x:x/geoserver/rest/"

headers1 = { "Content-type" : "application/xml", "Authorization": authorizationHeader}
coverageBaseURL = baseURL + "workspaces/Radar/coveragestores?configure=all"

headers2 = { "Content-type" : "image/tiff", "Authorization": authorizationHeader}
uploadBaseURL = baseURL + "workspaces/Radar/coveragestores/"

for file in os.listdir(resultPath):
if file.endswith(".tif"):
    name = os.path.splitext(file)[0].split("_")[4]

    #Create coverageStores for each file in GeoServer
    data1 = '<coverageStore><name>' + name + '</name><workspace>Radar</workspace><enabled>true</enabled><type>GeoTIFF</type><url>/radar/</url></coverageStore>'
    req1 = urllib2.Request(coverageBaseURL, data1, headers1)
    response1 = urllib2.urlopen(req1)
    the_page1 = response1.read()

    #Upload files to GeoServer
    url2 = uploadBaseURL + name + "/file.geotiff"
    data2 = open(file, 'rb').read()
    req2 = urllib2.Request(url2, data2, headers2)
    req2.get_method = lambda: 'PUT'
    response2 = urllib2.urlopen(req2)
    the_page2 = response2.read()

.. I get the following error in GeoServer's log:

04 Nov 14:06:02 ERROR [config.impl] - Error occurred processing a configuration
change listener
org.geoserver.catalog.CatalogException: Error persisting GeoServerInfoImpl[null]
 to global.xml
...
Caused by: java.io.IOException: Failed to move %GEOSERVER%\data\global.xml.3a26dbde-4e4f-4f35-8fba-be3f81e2debb.tmp -
unable to remove existing: %GEOSERVER%\data\global.xml

So the error must be related to the loop, somehow.


texts: Updating an existing field with field calculator in QGIS modeler

texts:
When incorporating the field calculator into the graphical modeler in QGIS is it possible to update an existing field without having to create a new field and then remove the existing.
For example I have a field called "MEAN" and lets say we have one record with a value of 100. Is it possible to update that field with the following simple sum "MEAN" / 10. I know this is possible outside of the modeler environment by using the "Update existing field" tick box but I can't seem to replicate this in the modeler


texts: Intersect / Select by Location with WMTS file

texts:
I have points and WMTS tiles from NOAA.
Is it possible to find the intersection of the points and the tiles or would I have to download the actual raster data from NOAA?
I ask because the rasters are obviously very large, and the tiles are not.
Is there any way to 'use' WMTS layers, or are they just for display?


texts:
I am facing a challenge with the encoding in my SLD file in GeoServer on a Ubuntu server.
When I create a SLD file containing the danish letters æ, ø, å in the GeoServer web interface like this filter:

<ogc:Filter>
 <ogc:PropertyIsEqualTo>
   <ogc:PropertyName>priority</ogc:PropertyName>
   <ogc:Literal>Høj</ogc:Literal>
 </ogc:PropertyIsEqualTo>
</ogc:Filter>

In the "Style Editor" the "ø" looks fine and when submit I can see (in Firebug) that the posted request contains the letter "ø".
But when I try to edit the style again, my filter looks like this:

<ogc:Filter>
 <ogc:PropertyIsEqualTo>
   <ogc:PropertyName>priority</ogc:PropertyName>
   <ogc:Literal>H?j</ogc:Literal>
 </ogc:PropertyIsEqualTo>
</ogc:Filter>

The "ø" has been replaced with a "?". My SLD header looks like this:

<?xml version="1.0" encoding="UTF-8"?>

I have tried to use Latin-1 encoding

<?xml version="1.0" encoding="ISO-8859-1"?>

But I the "ø" still shows as a "?".
When I run env | grep LANG in the terminal on the Ubuntu server I get this output:

LANG=en_GB.UTF-8
GDM_LANG=en_GB
LANGUAGE=en_GB:en

This should indicate the default encoding is utf-8.
When I run a file --mime-encoding myld.sld I get:

us-ascii

My question is:

How do I get GeoServer to write the sld file in utf-8 encoding or in Latin-1 so I can use the danish letters in my SLDs?


texts:
How can make topology at nodes, link?

if nodes layer a and links layer b,
they have correlation each other.

then how can make topology for them?

for example, edit some nodes' location then change links location automatically.


texts: How to convert PostGIS WKB to Shapefile - Fiona 'Polygon' object not iterable error

texts:
I am relatively new to Python and programming and trying to create a shapefile using multiple WKBs queried from a PostGIS database. I am using psycopg2 to query the database and using shapely and fiona for format conversion. The output I am looking for is a shapefile with multiple polygons (rows) based on how many WKBs were queried. However, I am getting an error: "TypeError: 'Polygon' object is not iterable".

import psycopg2
from shapely.wkb import loads
from shapely.geometry import MultiPolygon, mapping
import fiona
from fiona.crs import from_epsg

def main():
    try:
        print("Connecting to  database...\n")
        db_connection = psycopg2.connect(host = "server.com", dbname = "database", user = "postgres", password = "postgres")
        print("Connected to database\n")
    except:
        print("Unable to connect to database")

    db_cursor = db_connection.cursor()

    try:
        db_cursor.execute("SELECT feature FROM table LIMIT 10")
    except:
        print("Unable to query")

    def wkb_to_shapefile(db_wkb, db_wkb_list, schema):
        for row in db_wkb:
            db_wkb_list.append(row)
            image_wkb = db_wkb_list[0][10]
            load_wkb = loads(image_wkb, hex = True)
            print(load_wkb)
            geometry = list(load_wkb)
            print(geometry)
        with fiona.open("output/shapefile.shp", "w", crs = from_epsg(4326), driver = "ESRI Shapefile", schema = schema) as output:
            output.write({"properties": {"feature": ""},
                          "geometry": mapping(geometry)})

    db_rows_list = []

    db_rows = db_cursor.fetchall()

    shapefile_schema = {"geometry": "MultiPolygon", "properties": {"feature": ""}}

    wkb_to_shapefile(db_rows, db_rows_list, shapefile_schema)

    db_cursor.close()
    db_connection.close()

if __name__ == "__main__":
    main()


texts: Getting GeoJSON file to load in OpenLayers

texts:
I have validated the JSON and I'm certain it is clean. I'm loading it from a web server (using Node.js).  The Firefox error is "Uncaught SyntaxError: JSON.parse: unexpected character at line 1 column 1 of the JSON data". Chrome reports the unexpected character to be "<".
It appears that the payload from the call to JSON.parse is the application's index.html file rather than the data file, but I may be interpreting this wrong.
Am I missing something obvious?
I've stripped the code back to as simple as I can for now:
index.js
import 'ol/ol.css';
import Map from 'ol/Map';
import OSM from 'ol/source/OSM';
import GeoJSON from 'ol/format/GeoJSON';
import VectorLayer from 'ol/layer/Vector';
import VectorSource from 'ol/source/Vector';
import View from 'ol/View';
import {fromLonLat} from 'ol/proj';
import TileLayer from 'ol/layer/Tile';

//CONFIG
const initmapcenter = [-124.01901771596816, 54.018354087405754 ];
const dataurl = './data/SABoundarygeo.json';

// Polygons style
function polygonStyleFunction(feature, resolution) {
    return new Style({
      stroke: new Stroke({
        color: 'blue',
        width: 1,
        }),
        fill: new Fill({
            color: 'rgba(0, 0, 255, 0.1)',
        }),
    });
}
// Vector layer from a GeoJson file url
var vectorPolygons = new VectorLayer({
    source: new VectorSource({
      format: new GeoJSON(),
      url: dataurl
    })
    // , style: polygonStyleFunction
  });

const map = new Map({
  target: 'map',
  layers: [
    new TileLayer({
      source: new OSM()
    }),
    vectorPolygons
  ],
  view: new View({
    center: fromLonLat(initmapcenter),
    zoom: 5
  })
});

index.html
<!DOCTYPE html>

    <meta charset="utf-8">
    <title>Using Parcel with OpenLayers</title>


texts: Create QGIS action that shows point at coordinate on click

texts:
In QGIS 2.18 in my current point-layer I have related coordinates, which I would like to pop up when clicking on a point in the point layer. I try to use QGIS "Action" for this.

Simplified example of my point-Layer:

ID_____ | X_____ | Y_____ | Xrel__ | Yrel__

----------

A______ | 1.0___ | 1.0___ | 1.2___ | 1.3___
A______ | 1.0___ | 1.0___ | 1.5___ | 1.0___
A______ | 1.0___ | 1.0___ | 1.0___ | 1.1___
B______ | 2.0___ | 2.0___ | 2.2___ | 1.9___

So when clicking on point A (or nearby coordinate of point A) I would like to highlight the coordinates of Xrel and Yrel corresponding to A. Is there a way to do this using the Generic or Python action type?

Preferably when selecting another point or deselecting, the highlighted coordinates disappear.

I made a sample script that follows the steps I want to do. It is not in pygis style, but more a guideline of what I want to achieve:

# here some code to load all necessary python gis modules

# create a small buffer around the point where the user clicked
coordinateBuffer = buffer(clickedPointOnMap,0.0001)

# find all features within that buffer in the feature class (point layer)
pointsInBuffer = intersect(pointLayer,coordinateBuffer)

# statement to make sure to remove previous created points
    if isempty(pointsInBuffer)
        if isnotempty(pointsOnMap)
        remove(pointOnMap)
    else break
    else
    if isnotempty(pointsOnMap)
        remove(pointOnMap)
    else
        # extract coordinates in the attribute table of the point layer
        toHighlight = pointsInBuffer[:, (Xrel,Yrel)]

        # create a point and set layout of the toHighlight points layer
        pointsOnMap = make(toHighlight)
        pointsOnMap(fill_color_rgb(255,211,0))}


texts: How To Create "Boating Time" Areas -- ArcGIS

texts: I'm fairly new to GIS, and I would like to create what I'm calling a "boat time area". Similar to a drive time area, it would give an area of all the places a boat can be in a given amount of time, with a given starting point. I need to account for terrain that boats would have to navigate around (islands, coasts, etc). How would I achieve this? What tools should I use?


texts: Add more vertexes along GeoJSON linestring

texts:
How can I add more vertex
Some lines only have two vertices. The problem is on 3D globe it's not dense enough and does not clamp to earth

Does ogr2ogr or OGRINFO support this geospatial operation


texts: Converting coordinate systems of shapefiles using QGIS?

texts:
I have two shapefiles representing the same GPS track taken at different times.  However, they have different coordinate systems so they do not match, they are separated by about 354 meters E/W.  I would like to convert them to the same coordinate system so that they overlay each other.  Both were exported directly from Trimble Pathfinder Office.

Shape1 has the CRS +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +units=m +no_defs

Shape2 has the CRS +proj=longlat +datum=WGS84 +no_defs

in QGIS I saved Shape2 to the CRS of Shape1, but it did not change its position.  I also saved Shape1 to the CRS of Shape2, and it also did not change its position.  In both cases the two shapefiles remained separate by about 354 meters.


texts: Search or filter values in two-segment list in Google Earth Engine

texts:
i have a two-segment list;
var list =ee.List([[101,12155],[102,12157],[103,12158]])

i want to search or filter for 102 and get result (12157), for 103 and get result(12158)
How can i do it ?
i find a way (below code) to get result but is there a simpliest way?
var list = ee.List([[101,12155], [102,12157], [103,12158]])

// get ui.select label
function findObjectByValue(slist, key) {
for (var i = 0; i < slist.size().getInfo(); i++) {
var number = ee.List(slist.get(i)).get(0)
print(number)
if (number.getInfo() == key) {
  //print('AAA')
  return ee.List(slist.get(i)).get(1);
  }
}
return null;
}

print(findObjectByValue(list, 101)) // print 12155


texts: How to get the last drawn point for a Geometry in OpenLayers

texts:
I wonder if I can get the position of the last point for a geometry in OpenLayers?
For example, I have an OpenLayers.Geometry.LineString, how to get its last point?
Clarification in the comments:
In fact, what I want is the last drawed point for the polygon. Is it impossible?


texts:
Has anyone been using callout labelling with QGIS, either trunk or stable?

Used for clustered points. That means the point stays where it is, the label maybe moved if required and a leader (plain or pointed) shows the linkage between the label and the point.

How could I do it (you know what I mean, no coding....)?


texts:
In my job are currently looking for a ArcGIS free alternative, we try QGIS but ArcGIS reshape feature  is very important for us when working with curves. What we do basically is:

We've a line as marked in the red square bellow;

then we select two point of that line as show in the next picture;

and the final result is the geometry reshaped as show in this final picture;

Is there any feature or plugin in QGIS for accomplish this behavior?


texts: How to label points with gradient-colored data-defined rectangles in QGIS 3.18?

texts:
I am trying to symbolize a point layer with rectangles shaped and colored by gradient according to attribute values. I am trying to do this via label backgrounds rendered as "Geometry generator" marker symbols to prevent rectangles from overlapping and enable easy callout creation (and to enable gradient fill -- this is, to my knowledge, impossible with Ellipse markers).
Here is an overview of what I am trying to achieve. The white rectangles (these are label background symbols) are placed as I want, but I cannot color them properly. The colored rectangles are created with geometry generator, they are the desired symbols, but I cannot place them correctly.

I encountered a feature of geometry generator I cannot explain. Let us take the following table as an example point layer:
x   y
10  12
15  27
43  43
27  15
12  10

When I pass the code $geometry to geometry generator in the Symbology section, the points are placed where they are meant to be:

However, when I add labels and pass the same code to geometry generator in Labels > Background > Marker symbol, the points are placed well away from the point locations:

These labels change their locations as I drag the map canvas, they enclose the point layers when the points are near the upper left corner of canvas. The locations of text labels remains as usual.
I have two (three) questions:
1. Why are the geometry-generated label background symbols placed like that?
2. How can one avoid such improper placement?
(3.) Maybe there is another way to symbolize points with geometry generator and prevent them from overlapping that is less painful?


texts: How to change filter and style of vector layer after initialized?

texts:
I want to change filter and style like this

-user enter search condition ==> filter condition change
-user enter color of feature ==> style color of feature change

i know how to add filter and style for vector layers when initialization this vector, but
i dont know how to change style after that.


texts: Different symbology for each theme in QGIS

texts:
I want that a layer changes its symbology depending on the theme selected.
Is that possible?
Maybe if I could retrieve the name of the active theme I could add it as a part rule at the symbology.
If ThemeName == 'A'
...do this symbology
else
...do this one

There's a a way to get it?


texts: Updating raster pixel values where they overlap with another raster band

texts: I am trying to replace/update the raster pixel values of band '3' of a Land Cover Map to a new value of '24', but only where these pixels overlap with value '4' of another raster layer. How do I write this as an expression in the raster calculator? I am using QGIS.


texts: M-N relationship class in geodatabase with existing linking table?

texts:
I have a parcel database with an M-N relationship (Parcels to Deeds). I have already created and populated a linking table. We have also set up two relationship classes in the geodatabase for the two 1-M relationships (Parcels-Links, Deeds-Links).

Now we want to query that in a webmap using the Esri Javascript API, so when they click on a parcel it lists the deeds. We know we can do this with two nested relationship queries, but it's cumbersome. We were hoping that we could set up the geodatabase so we could query it in one shot. It seems common enough that Esri would have a system for it.

We first tried creating a view that joins the links to the deeds, so we could do one relationship class from parcels to the view, but it will not allow us to link a relationship class to a view (is that correct?).

Then we tried creating an M-N relationship class, but it insists on creating a new linking table and I can't find an option for using an existing table.

Any ideas?


texts: Python Code: Help with Printing a Message from an AddMessage Result?

texts:
I have written a code to count sheets or pages in an mxd. How would I print a message before the display of the count result? For Example: the count is 11.  I want it to display "Sheet Count : 11". My code is below:

import arcpy, os, sys, string
from os import path as p
from arcpy import mapping as m

Make parameters for people to choose mxd and folder of PDFs to compare

mxdList = arcpy.GetParameterAsText(0).split(";")
pdf_path = arcpy.GetParameterAsText(1)

Use Search Cursor to go through Attribuite Table to get Sheet number info

count = 0
mxd_dict = {}
for mapDoc in mxdList:
    arcpy.AddMessage(mapDoc)
    mxd = arcpy.mapping.MapDocument(mapDoc)
    max_list = []
    for lyr in m.ListLayers(mxd):
        try:
            rows = arcpy.SearchCursor(lyr)
            for row in rows:
                max_list.append(row.Sheet_ID)
                count += 1
            sheet_count = max(max_list)
            mxd_name = p.basename(mapDoc).split('.')[0]
            mxd_dict[mxd_name] = sheet_count

        except:
            pass
arcpy.AddMessage (count)
 Right here is probably where I need the code. Sheet Count :


texts: Splitting a branched river network at point layer

texts:
I need to split a country-scale river network (polyine) at certain barriers (point) so the length of (branched) river sections that have a certain minimum length. I used the following procedure:

1) Merge/Unsplit river network to one long line (only one row in the attribute table).

2) Split line at Point.

The problem is now that it not only splits the lines at the barriers, but also at intersections/branchings within the river network. I also tried Xtool's split line feature with similar results, and gvSIG's sextante toolbox didn't work at all for me. Any suggestions?

I am using ArcGIS 10 SP5, but also open to solutions in Qgis or gvSig


texts: How to create a symbol with no fill color in qgis?

texts: Where can I find a simple circle to draw it around a point feature. I can change symbol fill color as well as symbol border color. But I didn't find a way to set fill color to no color.


texts: How to remove all transparency from all layers?

texts: I am using ArcMap 10.1 and I would like to remove all transparency from all the layers in the map at once.


texts: Pseudo code to correct Geographic points for UTM Easting distortion on Web Mercator Projection

texts:
I have chosen to use a Web Mercator projection system for:

Use of Tiling system for low memory usage
Accessing only one site/area at a time
Reuse of Google Tiles - for many sites

In some cases I will use high resolution UTM images, projected onto the tiled Web Mercator user interface. I suspect it is due to Easting distortion increasing in severity at greater distances from the central meridian.

see Calculating areal distortion outside UTM zone?

This is inaccessible to me as a Software Engineer with inadequate experience in GIS. Needing:

To know if that linked answer has the full capability for my situation
To know how to represent that as two Pseudo code functions: i) Point GeoToUTMPoint(Point Geo); and ii) Point UTMToGeoPoint(Point UTM).

You need to download the images and toggle between them to see the distortion. The site here is in Geelong, Australia with UTM Zone 55S. Markers placed on the map (see second UTM Tiled image) and have an assumed Geo-Coordinate. When those are placed in Google Maps there is an offset error. The top left and bottom right are quite stable, areas to the top right and bottom left increasingly shift the coordinates generally anti-clockwise but moreso shearing than fully rotating.

I know there's a reason Web Mercator is resisted being non-conformal. But I'm hoping there are some here who can at least point me in the right direction with a few pseudo code functions that I can use together.

I'll pay for the answer and keep it published for all: https://www.upwork.com/jobs/Function-for-Correcting-Distortions-from-UTM-Projection_~01661dc8f1a19fbb8c. If you want the money you need to prove you answered it of course, which can happen if you publishing a code (I supply from upwork chat) in a comment with your same StackExchange alias.


texts: Can not Run function when item selected in comboBox. Pyqgis. QGIS

texts:
I'm writing a plugin, I want to change the content of a comboBox (fields) depending of the content of another comboBox (layers).

My plugin is a dockwidget. I can populate the fields comboBox just starting the plugin in the def run(self): function, and it's working fine:

        layerNom = self.dockwidget.comboBoxErrores.currentText()#Get the text value of the comboBox
        for lyr in QgsMapLayerRegistry.instance().mapLayers().values():    #Get layers
            if lyr.name()==str(layerNom):
                fields = lyr.pendingFields() #Get Fiels
                field_names = [field.name() for field in fields] #Field List
                self.dockwidget.comboCampos.addItems(field_names) #Added to the comboBox

But when I try to update the content of the fields comboBox I really don't know how to do it.

In my "pluginName_dockwidget.py" file, in the init function I have created the connection with the function that will repopulate the comboBox with the new content:

self.comboBoxErrores.currentIndexChanged['QString'].connect(self.pueblaCampos)

In the same file I have a function that get the value of the layer name of the comboBox and emit it with the signal created previously:

sennalComboIndex = pyqtSignal(str)
.
.
.
def pueblaCampos(self):

    layerName = self.comboBoxErrores.currentText()
    self.sennalComboIndex.emit(layerName)

In my main file, I have the rest of the function:

def pueblaCampos(self, layerName):

    for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyr.name()==str(layerName):
            fields = lyr.pendingFields()
            field_names = [field.name() for field in fields]
            self.dockwidget.comboCampos.clear()
            self.dockwidget.comboCampos.addItems(field_names)

When I change the slected item in the layers comboBox, nothing happens. What Am I missing?


texts: I am very much fine with creating a circle in Geoserver using st_buffer to create a view by providing latitude, longitude and radius values. It creates only a circle but without center. But I want to even see the center of the circle marked. What is to be done for this???