texts: List layers in dataframe as parameters in script validation

texts:
I'm trying to do some script validation in a tool I wrote. The user first selects a dataframe from the current mxd. This parameter works fine. When the user makes this selection, the next parameter, a list of layers in that dataframe should populate the next parameter. But all I get with the code below is an empty list. I can get layers for the entire mxd, but am unable to restrict it to the selected dataframe. This seems like something simple I'm missing. How can I accomplish this?

def initializeParameters(self):
"""Refine the properties of a tool's parameters.  This method is
called when the tool is opened.
Populates data frame parameter with data frames in map and disables other
parameters.
"""
try:   # this works
    mxd = arcpy.mapping.MapDocument('CURRENT')
    dflist = [df.name for df in arcpy.mapping.ListDataFrames(mxd)]
    self.params[0].filter.type = "ValueList"
    self.params[0].filter.list = dflist
    del mxd
except:
    arcpy.AddError("\nMust run this tool in ArcMap, not ArcCatalog")

self.params[1].enabled = False
self.params[2].enabled = False

return

def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""

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

    if self.params[0].value:
        self.params[1].enabled = True

        # this does not work the way I want
        layerList = [lyr.name for lyr in arcpy.mapping.ListLayers(mxd, self.params[0].value) if not lyr.isGroupLayer]

        uniqueList = list(set(layerList))
        uniqueList.sort()

        self.params[1].filter.type = "ValueList"
        self.params[1].filter.list = uniqueList


texts: Fill shapefile column with raster value, shapefile with alias name and calculated field

texts:
I've to sample raster values to fill a column of point shapefile using Python, and for this I've no problem. But this shapefile has some fields with alias name and some fields that result from an expression.
How can I copy this fields to a new shapefile and then fill the column with raster values?
I have created the new fields and using SetField from ogr to copy the values of original shapefile but I can't found a way to set the Alias name and the calculated fields.
This is my code, I'm not a programming expert:
def sample(shp,raster):
#open raster
        rasterfile = raster
        raster = gdal.Open(rasterfile)
        gt = raster.GetGeoTransform()
        banda = raster.GetRasterBand(1)
#open shapfeile
        driver = ogr.GetDriverByName('ESRI Shapefile')
        dataset = driver.Open(shp)
        layer1 = dataset.GetLayer()
        spatialRef = layer1.GetSpatialRef()
#create new shapefile
        outputShapefile = shp.split('.')[1] + '_coord_z.shp'
        print(outputShapefile)
        driver = ogr.GetDriverByName('ESRI Shapefile')
        new_data = driver.CreateDataSource(outputShapefile)
        layer = new_data.CreateLayer(outputShapefile, spatialRef,ogr.wkbPoint)
#define fields name
        name_field = ogr.FieldDefn('Apparati', ogr.OFTString)
        name_field.SetWidth(100)
        layer.CreateField(name_field)

        name_field = ogr.FieldDefn('Stato', ogr.OFTString)
        name_field.SetWidth(15)
        layer.CreateField(name_field)

        name_field = ogr.FieldDefn('coord_Z', ogr.OFTReal)
        layer.CreateField(name_field)

# sampling raster values
        for feat in layer1:
                geom=feat.GetGeometryRef()
                point=geom.Centroid().ExportToWkt()
                mx=float(point.split(' ')[1][1:len(point.split(' ')[1])])
                my=float(point.split(' ')[2][0:len(point.split(' ')[2])-1])
                px = int((mx-gt[0])/gt[1])
                py = int((my-gt[3])/gt[5])
                intval = banda.ReadAsArray(px,py,1,1)
                zcoord=float(intval)
                #set layer field's values
                feature = ogr.Feature(layer.GetLayerDefn())
                feature.SetField('Apparati',feat.GetField('Apparati'))
                feature.SetField('Stato',feat.GetField('Stato'))
                feature.SetField('coord_Z',zcoord)
                feature.SetGeometry(geom)
                layer.CreateFeature(feature)
                feature = None
        raster = None
        Dataset = None
        new_data = None

I need to add to fields 'Apparati', 'Stato', and 'Coord_Z' an alias name that can be show in attribute table, and for calculated fields I've no idea of what I can do. In this example I show just 3 fields but I have much more in my shapefile. Do you think it is a good idea to use OGR/GDAL or there are some better libraries?
P.S. The raster I need to sample is a DEM
P.P.S. Even if I use the QGIS processing tool "Sample raster values" I've the same problem with alias names and calculated fields


texts: Possible to execute external python script from ArcMap Python Window?

texts:
I have this python script that I execute inside the ArcMap Python Window and it does an array of spatial analysis processes to get me to my desired final mxd content. I then have another script that overrides my existing feature classes on my AGOL account with the current layers in the open mxd. : Click Here (it is necessary that this script is ran outside of ArcMap via the Python shell. )

The only issue is that after I run my first script in ArcMap I then have to open my python shell outside of ArcMap in order to run the second script. I am trying to find a way to execute the 2nd script from the first script without manually leaving arcmap to open my python shell and load the script. But keep in mind my first script is being run from the ArcMap Python Window. SO I cannot figure out a way to escape that... hope this question makes sense and is this even possible. I am basically trying to eliminate all user-interaction to completely automate this task once I load the first script into ArcMap.

Please feel free to ask if this does not make sense.


texts:
I have a GeoTiff that loads incorrectly in QGIS.

Instead of displaying the image at the correct coordinates, QGIS positions the raster with its upper left corner at (0,0).

This same file loads correctly in Google Earth, but for some reason QGIS isn't recognizing the georeferencing info.

The original motive was to load USGS topo quadrants into QGIS.  I downloaded the georeferenced PDF files available from USGS, then used GDAL to convert:

gdal_translate -a_srs EPSG:32617 example_USGS_topo.pdf example_USGS_topo.tif --config GDAL_PDF_LAYERS "ALL" --config GDAL_PDF_LAYERS_OFF "Images" --config GDAL_PDF_DPI 300 --config GDAL_PDF_BANDS 3

More info on that in the helpful USGS instructions here: Converting US Topo GeoPDF Layers to GeoTIFF

If you'd like to work with the same USGS PDF file I'm using, here's a link to the original USGS ZIP containing the PDF. https://ims.er.usgs.gov/gda_services/download?item_id=8139809

Note that the above gdal_translate command applied to this example PDF produces a GeoTiff with the following gdalinfo output:

Driver: GTiff/GeoTIFF
Files: NC_Ashford_20160624_TM_geo.tif
Size is 6825, 8700
Coordinate System is:
PROJCS["WGS 84 / UTM zone 17N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-81],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32617"]]
GeoTransform =
  408447.1636098031, 2.03190732774024, -0.01942631127630456
  3972082.073249529, -0.01942631127630456, -2.03190732774024
Metadata:
  AREA_OR_POINT=Area
  AUTHOR=USGS National Geospatial Technical Operations Center
  CREATION_DATE=D:20160624090308Z
  CREATOR=ESRI ArcSOC 10.0.2.3200
  KEYWORDS=Topographic, Transportation, Hydrography, Orthoimage, U.S. National G
rid, imageryBaseMapsEarthCover, Imagery and Base Maps, Geographic Names Informat
ion System
  NEATLINE=POLYGON ((421074.64910106914 3970498.3348052166,420940.92899800668 39
56511.7959603746,409525.75346028153 3956620.932213902,409659.47356334399 3970607
.4710587487,421074.64910106914 3970498.3348052166))
  SUBJECT=This image map depicts geographic features on the surface of the earth
.  It was created to provide a representation of accessible geospatial data whic
h is readily available to enhance the capability of Federal, State, and local em
ergency responders for homeland security efforts.  This image map is generated f
rom selected National Map data holdings and other cartographic data.
  TITLE=USGS 7.5-minute image map for Ashford, North Carolina
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  408447.164, 3972082.073) ( 82d 0'51.74"W, 35d53'19.42"N)
Lower Left  (  408278.155, 3954404.479) ( 82d 0'51.18"W, 35d43'45.67"N)
Upper Right (  422314.931, 3971949.489) ( 81d51'38.62"W, 35d53'19.44"N)
Lower Right (  422145.922, 3954271.895) ( 81d51'39.16"W, 35d43'45.67"N)
Center      (  415296.543, 3963176.984) ( 81d56'15.18"W, 35d48'32.64"N)
Band 1 Block=6825x1 Type=Byte, ColorInterp=Red
Band 2 Block=6825x1 Type=Byte, ColorInterp=Green
Band 3 Block=6825x1 Type=Byte, ColorInterp=Blue


texts: OSRM 4.8 import planet-latest.osm.pbf 31G

texts:
I use OSRM version 4.8.

My VM configuration:

16 CPU
64 RAM
150G SWAP
150G stxxl disk size
250G data directory
50G root

When I import Europe-latest (17Gb) file - it takes about 6 hours to start routing.

When I import Planet-latest (31Gb) file - it takes more than 4 days! And it is still on the preparing step. All disk activities are only on swap disk.

When I start extract step it uses all 16 CPUs. But the preparing step uses only 1 (ignoring '-t 16' option)

Is this the expected performance?


texts: Fill no data pixels for specific area (not the whole image) resulting from cloud masking of Landsat 8 multi temporal composite (Google Earth Engine)

texts:
I need to conduct a Landuse classification of New Caledonia. I therefore used the Landsat 8 satelite images of one year and created a multi temporal mosaic by a mean function.
As most of the areas were ostly clouded, I used a cloud mask, which resulted in areas of no value.
I then tried the focal.median function. However, this led to a change of DN values for the whole image.
I am therefore looking for a function which only tackles the specific areas of no data and leaves the other values untouched.
I was thinking of using Sentinel-1 Radar data, check the radar signature of those areas of my landsat8 image where I have no values, look for the same radar signature at another area, from where I also have a spectral value in my Landsat8 image to then use that value of my Landsat 8 image to fill my missing Landsat8 areas.
Does that make sense? If so, how could I execute that idea? Or is there even any other option?


texts: On arcmap 10.1, I have elevation data that is shown as contour lines (polylines) i also have point data showing the distribution on certain animals. I want to run some kind of statistical analysis to see at what height the animals are most commonly found at, how would i do that/what tools would i use? Thanks


texts: Publishing individual services without ArcGIS for Server?

texts:
I have many layers in my map that I want to publish as individual services. Meaning they all have their own rest end point, and are not grouped into the same hosted feature service. I cannot use ArcGIS Server to complete this task.  Each item will need it's own metadata and tags in the future.

What is the best way to publish these items? Some are group datasets, too.  I would publish those as a single hosted feature service so they remain under the same metadata and tags and rest end point.

Should I individually publish them through the "Publish a Service" tool? Should I publish a map package and then publish the data through online? Should I create layer packages? Should I create and ArcPy script like this -> Overwrite existing map service in ArcGIS Server 10.1 (through ArcPy) ? What is the best way to increase efficiency, conserve credits, and ensure that the data is transferred properly?

This is more of a discussion than a true question, but I want to know the opinions of the community? Many of my clients cannot afford server and we work strictly with cloud-hosted services. Perhaps there is a tool I am not aware of to help increase efficiency when publishing large amounts of data?


texts: MapBox GL Native SDK iOS - SnapShotter - remove watermark

texts:
https://docs.mapbox.com/ios/maps/examples/map-snapshotter/
https://docs.mapbox.com/archive/ios/maps/api/5.3.0/Classes/MGLMapSnapshotter.html
Trying to figure out if there is a way to not include mapbox watermark in MapBox GL Native SDK for iOS


texts: How to use Replace function in a custom script?

texts:
I would like to use Replace function in a QGIS model but it doesn't work as I wish. I was working with Joseph's brilliant solution (using a custom script in a model) but I couldn't find out how to replace a specific character in a field. I would like to replace every '\' in Fromfield with '/' and write the result into Tofield.
I tried this and a few more, but nothing:

##Test=name
##Layer=vector
##Tofield=Field Layer
##Fromfield=Field Layer

from qgis.core import QgsExpression

layer = processing.getObject(Layer)
layer.startEditing()

idx = layer.fieldNameIndex(Tofield)
g= str(Fromfield)
e = QgsExpression(g.replace("\\","/"))
e.prepare(layer.pendingFields())

for f in layer.getFeatures():
    f[idx] = e.evaluate(f)
    layer.updateFeature(f)

layer.commitChanges()

Is this possible at all?


texts: Correcting Landsat images shifted by strangely high distance using ArcMap?

texts:
I have a set of ~20 Landsat images acquired with sensor version varying from LT5 up to L8, downloaded from EarthExplorer. Many of them do not align to the most current images, they are shifted by a distance between 100m (what is understantable) to 500m or a few kilometers. But the older ones have the misplacement even bigger - one of the LT5 images is shifted by 40 km! (though the corelation between error scale and Landsat version is none for me, I also have LT5 images with a very small displacement)

Using ArcMap, how can I fix this to have all images in the same place? I need to be 15m precise in my analysis (pansharpened 30m).


texts: QGIS Server Internal Server Error getcapabilities request

texts:
I followed the links from http://anitagraser.com/2012/04/06/qgis-server-on-windows7-step-by-step/ for installing the QGIS server at my Win 7 64bit system but the getcapabilities request doesn't work.

That's why I try to start the qgis_mapserv.fcgi.exe directly. It's said that many libraries are not found. I added all libraries from osgeo/bin to /osgeo/apps/qgis/bin.

Now the qgis_mapserv.fcgi.exe starts and give the output:

  QFSFileEngine::open: No file name specified

The apache error-log doesn't show anything.
The getcapabilities request doesn't work and
I get the following error-text in the browser:

  Internal Server Error

The server encountered an internal error or misconfiguration and was unable to complete your request.

  Please contact the server administrator, admin@local.domain and inform
  them of the time the error occurred, and anything you might have done
  that may have caused the error.

  More information about this error may be available in the server error log.

I read some questions with the same problem here but in my case they doesn't work.

I added all librarys
installed the mod_fcid separately from osgeo4w

Has somebody any other ideas?


texts: How to get estimates of areas of overlapping polygons in one shp to another?

texts:
Sorry for my ineptitude, but I have two vector shape files of Africa. One has polygons of countries, and one has polygons of tribes a long time ago. What I want to do is to find the intersection of these two polygons, and get some idea of how many tribes were contained in a country today, and what share of land in the country today that they took up. Hopefully that makes sense. So, I'd like to end up with some data with a row for a country, say, Nigeria, that has:

NIGERIA    | Tribe 1: 20% of land  |  Tribe 2: 40% of land  | Tribe 3: 40% of land |

Something along those lines. Is this possible in QGIS or ArcGIS? My preferred software is QGIS.


texts: I have a corrupt Shapefile. But, All other files are still good. How can I create or generate a Shapefile from only the ".shx .prj .mnd .drw .dbf"?


texts: Projects using compatible Sentinel-2 and Landsat-8 reflectance data

texts:
Right from the beginning on the Senintel-2 Satellites were constructed to harmonize with the Landsat-8 products in terms of their spatial and spectral capabilities. What I wonder is, if there are projects which really used this advantage. I know there is the NASA HLS Product. But I don't know of many projects that used the data produced by this campaign.
If anyone knows of more projects, I'd be curious to hear about them


texts:
Ho do I go about getting GDAL 2.0 on Windows? I can get the source from the OSGEO site http://www.osgeo.org/node/1591 but I have no idea how to compile it.

Is there not a built version of GDAL 2.0 like there is for the 1.X releases? I couldn't find an install option in the Osgeo4w installer either.


texts: How can I analyse the visibility of the urban area in ArcMap?

texts: I want to analyse the visibility of the urban area using ArcMap. I saw the image below and I want to have the same output. Can you explain the process and steps to obtain this output?


texts: WKID issue when adding raster layer from ArcGIS "mapserver" to QGIS via REST

texts:
QGIS docs here describe how to do this. I am trying to add a layer from this MapServer:

gps.digimap.gg/arcgis/rest/services/JerseyBaseMap/MapServer

When I try the following command from the Python console I get the error below.

>>> qgis.utils.iface.addRasterLayer("http://gps.digimap.gg/arcgis/rest/services/JerseyBaseMap/MapServer?f=json&pretty=true","raster")

If I instead try the gdal_translate command line option it appears to be a issue with the the Well Known Text ID not being recognised:

$ gdal_translate --debug on "http://gps.digimap.gg/arcgis/rest/services/JerseyBaseMap/MapServer?f=json&pretty=true" JerseyBaseMapLayer.xml -of WMS
HTTP: Fetch(http://gps.digimap.gg/arcgis/rest/services/JerseyBaseMap/MapServer?f=json&pretty=true)
WMS: Did not get WKID
HTTP: Fetch(http://gps.digimap.gg/arcgis/rest/services/JerseyBaseMap/MapServer?f=json&pretty=true)
ERROR 4: `/tmp/file.dat' not recognised as a supported file format.

GDALOpen failed - 4
`/tmp/file.dat' not recognised as a supported file format.

Question: As I know this provider's SRS (EPSG:3109) can I 'tell' QGIS/GDAL to apply it to the source data file, say by pointing to a WKT definition file?

gdal_translate has an option [-a_srs srs_def] to overide the output file projection, but I can't see a way to apply an SRS to an file imported via the REST API..

Very frustrating, as the data can easily be viewed in the ArcGIS webmap viewer:

www.arcgis.com/home/webmap/viewer.html?basemapUrl=http%3A%2F%2Fdcw7.digimap.gg%2Farcgis%2Frest%2Fservices%2FJBM%2FMapServer&source=sd


texts: Finding route between 2 points based on same attributes using ArcGIS Network Analyst?

texts:
I'm using Arcmap 10.5 and I would like to solve a network analysis problem.

I have a road network layer and two point layers. The first points layer contains the start points, the second one contains the end points (the same count). Both layer has route_id field. I would like to find routes based on this common route_id value.
So, I would like to find routes from each start point to each end point where the route_id values are the same (overall 100 routes).
They would be very-very simple routes, without any restrictions, barriers etc. Only the shortest paths are needed.

I tried this:

I've made a Network Dataset and loaded the road network layer.
After that I created a Closest Facility Layer and loaded start points as locations.
Next, I loaded end points as Incidents.

I've made routes successfully, but I would like to make routes based on the same attribute values, for example: from the start point with route_id=1 to the end point with route_id=1 etc. with 1 to 1 join.

How can I solve this problem?


texts: Plotting WFS layer from ArcGIS Online with OpenLayers?

texts:
I would like to plot WFS layer with OpenLayers.

I tried with a sample WFS stream at first and succedeed.

WFS stream : https://ahocevar.com/geoserver/wfs?service=WFS&version=1.1.0&request=GetFeature&typename=osm:water_areas&outputFormat=application/json&srsname=EPSG:3857

So now I want to adapt it my my FWS stream.

My WFS stream is from Atmo Hauts de France :
https://dservices8.arcgis.com/rxZzohbySMKHTNcy/arcgis/services/alrt3j_hdf/WFSServer?service=wfs&request=getcapabilities
I adapted it with getFeature instead of getCapabilities to get the features I want :
https://dservices8.arcgis.com/rxZzohbySMKHTNcy/arcgis/services/alrt3j_hdf/WFSServer?service=wfs&request=getFeature&typename=alrt3j_hdf:alrt3j_hdf_0
For the example above the output format is geojson, but unfortunately ArcGIS does not seem to produce WFS with geojson output format.

So I tried to use the GML format output. I don't get error message anymore ! But still nothing on the map ...

Here is the code sample that works with the WFS I used as example :

//map
  fond = new ol.layer.Tile({
        source: new ol.source.OSM()
      })

// WFS stream
  var waterNorthAmerica = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(extent) {
      return 'https://ahocevar.com/geoserver/wfs?service=WFS&' +
          'version=1.1.0&request=GetFeature&typename=osm:water_areas&' +
          'outputFormat=application/json&srsname=EPSG:3857&' +
          'bbox=' + extent.join(',') + ',EPSG:3857';
    },
    strategy: ol.bbox
  });

  var style_ligne_eau = new ol.style.Style({
      stroke: new ol.style.Stroke({
        color: 'rgba(0, 0, 255, 1.0)',
        width: 2
      })
    });
          // couche vecteur à partir du flux WFS
  var vectorLayer = new ol.layer.Vector({
    source: waterNorthAmerica,
    style: style_ligne_eau
  });
  //carte
      var map = new ol.Map({
    target: 'map',
    layers: [
      fond,vectorLayer
    ],
    view: new ol.View({
      center: ol.proj.fromLonLat([37.41, 8.82]),
      zoom: 4
    })
  });

Here is the code sample that plot nothing :

//map
  fond = new ol.layer.Tile({
        source: new ol.source.OSM()
      })
//wfs stream
    var url_hdf_test = "https://dservices8.arcgis.com/rxZzohbySMKHTNcy/arcgis/services/alrt3j_hdf/WFSServer?service=wfs&request=getFeature&typename=alrt3j_hdf:alrt3j_hdf_0&outputFormat=gml2&srsname=EPSG:3857";

    var hdf = new ol.source.Vector({
    format:  new ol.format.GML({
            srsName: 'EPSG:3857'
        }),
    url: url_hdf_test,
    strategy: ol.bbox
  });

  var style_point = new ol.style.Style({
      stroke: new ol.style.Stroke({
        color: 'rgba(0, 0, 255, 1.0)',
        width: 2
      })
    });
       var vectorHDF = new ol.layer.Vector({
    source: hdf,
    style: style_point
  });

  //plot the layers
      var map = new ol.Map({
    target: 'map',
    layers: [
      fond,vectorHDF
    ],
    view: new ol.View({
      center: ol.proj.fromLonLat([37.41, 8.82]),
      zoom: 4
    })
  });

For me both are similar, except the output format. If something see why the second one plot nothing that would be great !

UPDATE WITH SOLUTIONS :
I edit the post just to clarify the solution as TomazicM & Mike both gave me one.

Solution 1 : Keeping the GML format, but using ol.format.GML2 and not ol.format.GML.

Solution 2 : Using GeoJSON format, but setting the EPSG when defining GeoJSON format: ol.format.GeoJSON({dataProjection: 'EPSG:3857'})


texts: SLD version="1.1.0", GeoServer 2.12.3 and graphic stroke spacing

texts:
I am trying to get the spacing or gap between arrows in GeoServer to work. The rendering in QGIS is perfect but in the GeoServer the parameters are not read/taken into account.
I tried to get around it using Drawing arrow on line using SLD of GeoServer? and Specifying distance between arrows in line SLD in GeoServer? as well as the cookbook about dash array and graphic stroke. Everything is fine in QGIS but in GeoServer the gap disappears.

This is my code for this particular filter that needs the arrows:

<se:Rule>
  <se:Name></se:Name>
      <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
        <ogc:Or>
          <ogc:PropertyIsLike wildCard="%" escape="\" singleChar="_">
            <ogc:PropertyName>highway</ogc:PropertyName>
            <ogc:Literal>motorway</ogc:Literal>
          </ogc:PropertyIsLike>
          <ogc:And>
            <ogc:PropertyIsLike wildCard="%" escape="\" singleChar="_">
              <ogc:PropertyName>highway</ogc:PropertyName>
              <ogc:Literal>trunk</ogc:Literal>
            </ogc:PropertyIsLike>
            <ogc:PropertyIsLike wildCard="%" escape="\" singleChar="_">
              <ogc:PropertyName>oneway</ogc:PropertyName>
              <ogc:Literal>yes</ogc:Literal>
            </ogc:PropertyIsLike>
          </ogc:And>
        </ogc:Or>
      </ogc:Filter>
      <se:LineSymbolizer>
        <se:Stroke>
          <se:GraphicStroke>
            <se:Graphic>
              <se:Mark>
                <se:WellKnownName>filled_arrowhead</se:WellKnownName>
                <se:Fill>
                  <se:SvgParameter name="fill">#ff7f00</se:SvgParameter>
                </se:Fill>
                <se:Stroke>
                  <se:SvgParameter name="stroke">#000000</se:SvgParameter>
                  <se:SvgParameter name="stroke-width">0.5</se:SvgParameter>
                </se:Stroke>
              </se:Mark>
              <se:Size>14</se:Size>
            </se:Graphic>
            <se:Gap>
              <ogc:Literal>250</ogc:Literal>
            </se:Gap>
          </se:GraphicStroke>
        </se:Stroke>
      </se:LineSymbolizer>
    </se:Rule>

I am wondering if I am missing something between QGIS and my export to SLD, or if Geoserver does not read this for some reasons. I found this on github but as it is from 2016 I was wondering if it's still part of geoserver limitations.


texts: How to find out the co-ordinate system using .prj file

texts:
I'm very new to GIS concepts. Is there a way or procedure of finding out whether the co-ordinate system is in UTM, latitude/longitude or in stateplane using .prj file?

My .prj file contains the following text:

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
  84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]


texts:
I have a raster data in the database in separate rows. I would like to select 2 rows and do a bit operation with them (like |, &) and store the result in memory as a query output. Is this possible?

My databse looks like this:

Data | id
raster1 | 0
raster2 | 1

and I would like to do something like:

SELECT raster1 & raster2 FROM data

and query result are all pixels from both rasters.


texts:
I am trying to add a FEMA Hazard WMS layer to QGIS. I am using QGIS version 3.4.11-Madeira. I have used a shapefile in the past, but the hope with the WMS is that it would auto-update.

I went to https://hazards.fema.gov/femaportal/wps/portal/NFHLWMS and copied the url https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHLWMS/MapServer/WMSServer?request=GetCapabilities&service=WMS

Then I clicked add new layer, WMS, new, I pasted the URL WMS url above and clicked connect. I get the following response.

  The server you are trying to connect to does not seem to be a WMS server. Please check the URL.

  Instead of the capabilities string that was expected, the following response has been received:Could not get WMS capabilities: error occurred while parsing reference at line 86 column 118

  This is probably due to an incorrect WMS Server URL.

I cannot figure out what I am doing wrong.


texts: how to modify every point in Postgis geometry?

texts: I need to do some transform on Postgis geometry, saying turn wgs84 into gcj02, I can use ST_DumpPoints to get all the points, but after modifying them, it seems not an easy and elegant way to join them back. I tried ST_Collect, it's not taking the path from ST_DumpPoints.


texts: How to use user made png/jpeg icons in QGIS 2.0?

texts: I can only add svg icons in qgis 2.0. How can i use png or jpeg or gif images as icons in qgis2.0?


texts: Merge polygons in the same shapefile based on attributes, ArcGIS or ArcPy

texts:
I have a shapefile of polygons. I do not want any polygons with an area less than 25. For polygons with an area less than 25 I want to merge them with another polygon with an area greater than 25 that also has the same values for att1, att2 and att3.

I was able to write some code which found the polygons with areas less than 25 but I am stuck after that. can figure this out in ArcMap or ArcPy.

I want to automate this, I understand it can be done one by one but.


texts: Disable layer checkbox if no geometry on that zoomlevel (source=wms)?

texts: I have a map application using GXP which displays different data coming from different WMS sources. Some geometry is only available on certain zoomlevels (available via MaxScaleDenominator on the GetCapabilities). It should be only possible to select an overlay layer when you are at that minimum zoom level, otherwise it should be disabled. Any idea on how to do this?


texts: ArcGIS- Band collection statistics tool correlation less than -1?

texts:
I am trying to work out the correlation coefficient using band collection statistics in ArcGIS 10 between 2 rasters (a sea surface temperature raster and a storm genesis density raster created using point density tool).

However it is giving me a correlation coefficient of less than -1 which is not possible (as stated on the band collection tool help page)?

Any ideas on why I'm getting this value?


texts: Trying to use Python code and user defined inline parameters to build a SQL expression

texts:
I am trying to build a tool where the end user can run a model and make a selection from a feature class based on 3 different fields, calculate the total shape length, & export the results to a csv file.  I have everything else figured out except for one part. I want the user be able to make a selection based on all three fields or based on one , two, or all three fields without them having to edit the selection query.

As you probably guess when I don't populate one of the user defined parameters it pops up with an error. I think the way to do this would be to use the "calculate value"in model builder and use Python to build a sql expression and use that expression for the selection.  I was thinking of having 3 user defined parameters with drop down list of field values and using inline variable substitution in the calculate value tool to create a sql expression based on what the user defined in the parameters.

This is the expression I have right now for them to select using all three field values:

"DISTRICT" = '%District%' AND "SURFACE" = '%Surface%' AND "JURISD" = '%Jurisdiction%'

I have a basic understanding of Python.  Any help or suggestions would be greatly appreciated. Thanks!


texts: List of small lists in Python: How to retrieve items from small list in a loop?

texts:
I would like to automate the process of making some certain feature layers.

I made a list of lists that goes like that: [[feature class, subtype],[feature class,subtype],...]

I'd like to loop through the list and make feature layers using the items in the smaller lists for output paths and queries. So far I wasn't able to retrieve them correctly.

bigNeighborList=[["ExtractA",0],["ExtractA",1], ["BuiltupA",0],["CoastA",3],["SwampA",0],["WatrcrsA",0],["WatrcrsA",1],["LakeresA",0]] #["LakeresA",2],["LakeresA",1],["WatrcrsA",2],["GroundA",0],["Landfrm1A",1],["CropA",1],["CropA",2],["GrassA",0],["TreesA",1],["TreesA",0], ["FirebrkA",0],["SwampA",1],["SwampA",2],["CropA",3]]
for itemA in bigNeighborList:
    bigFC = [itemA[0] for itemA in bigNeighborList] #retrieve the first item - feature class name from a smaller list
    subtype = [itemA[1] for itemA in bigNeighborList] #retrieve the second item - subtype number
    outputItemA = os.path.join(wrksp100, str(bigFC[0]))
    print(outputItemA)
    qry = "%s = %d " % (subField,int(subtype[0])) ##query to choose the right subtype
    tempLyr = arcpy.MakeFeatureLayer_management(outputItemA,str(bigFC)+"_lyr",qry) #temporary layer for every list item

the print statement returns the path "C:\DAR...3\ExtractA" 8 times for every itemA in the list and creates new ExtractA tempLyr  (though named with all first items from small lists - "['ExtractA', 'ExtractA', 'BuiltupA', 'CoastA', 'SwampA', 'WatrcrsA', 'WatrcrsA', 'LakeresA']_lyr") over and over.

Is something wrong with my indexing?

I'm stuck.


texts: Get minimum distance direction from point to line

texts:
Is there a method using QGis to calculate the direction (degree) of the minimum distance from a point to a line.

I already have:
- Point shapefile with hundreds of points;
- Linear Shapefile with coast line;

I'm finding:
- New shapefile with minimum direction from to coastline with a column in its attribute table showing the degree (from north) to these lines.


texts: GEE create imagecollection of a specific tile

texts:
I'm trying to create an image collection from sentinel 2 chich will display aimages between two dates from a specific tile (T48QTD)

I couldn't find how to choose one specific tile .

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-10-14', '2019-12-08')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};


texts: Converting Polyline to Point Shapefile using arcgisscripting?

texts:
I would like to convert the vertices of an existing Polyline-Shapefile into a Point Shapefile, but I have a problem with reading out the geometries of the polylines (getPart()).

    outpath = os.path.split(SHPInput)
    Temp_PointShp = gp.CreateFeatureclass(outpath[0],"SohlPoint_temp.shp", "POINT",None, "ENABLED", "ENABLED") # empty PointShapefile
    gp.addfield_management(Temp_PointShp, "Z_Value", "DOUBLE")
    gp.addfield_management(Temp_PointShp, "M_Value", "DOUBLE")

    Srows = gp.SearchCursor(SHPInput)   # Cursor
    Irows = gp.InsertCursor(Temp_PointShp)

    row = Srows.next()
    while row: # for every row of the Polyline-Shapefile
        Geom = row.shape # Geometry
        OArray = Geom.getPart() # Array with Array

        PArray = OArray.getPart(0) # Array with Points

        Points = OArray.next()
        while Points: # for all vertices...

            row = Irows.NewRow() # create new row in Point-Shapefile

            row.X = Points.X
            row.Y = Points.Y
            row.Z = Points.Z
            row.M = Points.M
            Irows.InsertRow(row)

            Points = PArray.next()

        row = Srows.next()


texts: QGIS > Concave Hull plugin > Shared Nearest Neighbor Clustering

texts:
I have been using the Concave Hull plugin in QGIS & the Shared Nearest Neighbor Clustering feature.

I am having issues with this plugin, each time I choose the option of "Number of Neighbours" it doesn't seem to work.

This is the same issue DrottTom explains in his post here,

Grouping points in QGIS?


texts: Missing PostgreSQL driver in rgdal: a long and sad story

texts:
This issue has been discussed few times and any of the read posts have not been the ultimate answer for my case:

I'm trying to load data from a PostGIS database into the R Console. While many packages in R are suitable for it, it seems to be that only rgdal is able to implement directly spatial capabilities.

The bad new is (as mentioned in many posts) the installation of rgdal from CRAN in windows does not include the driver PostgreSQL (the driver is not listed by the command ogrDrivers() in the R Console).

I loaded the source for rgdal from the CRAN site to build the package myself as recommended in some posts but I always get an error message:

OGR_write.cpp:1:24: fatal error: cpl_string.h: No such file or directory

The same error appears when trying to install the package from GitHub:

library(devtools)
install_github("cran/rgdal")

Is there a way or a guide (for dummies) that can help me getting the PostgreSQL driver enabled in rgdal?

Further information that may be required to answer this question:

My operative System: Windows 7 64-bit
R-3.3.3 installed at C:\Program Files\R\R-3.3.3
GDAL-2.0.1 installed at C:\Program Files\GDAL
FWTools-2.2.8 installed at C:\Program Files (x86)\FWTools2.2.8
Rtools-34 installed at C:\Rtools


texts: How To calculate raster land cover inside overlapping buffer zones?

texts:
I am working with a raster land cover file that consists of several  land cover types (e.g. forest, buildings, grass ect.) in total there is 40 different land cover types.

I have 1200 bufferzonen (vector file) each on 500 meters. I would like, for each bufferzone (they have all a unik_ID) to calculate how many square meters that consist land cover type 1 (forst), 2 (buildings), 3 (grass) ect.

Some buffer zones are overlapping....

In the end I would like one crosstable that contains all the 1200 bufferzones and the square meters of each landcover type inside all the buffers.

I´m using ModelBuilder.


texts: Error when converting OSM relation type to GeoJson using osm2geojson

texts:
I have a function which is able to retrieve OSM data with Overpass API (Json) and convert the Json output to Geojson.
It works perfectly with ways, nodes and routes but it doesn't work with relation type.
Here my code:
def lines_OSM(query, lat, lon, radious, output, name):
overpass_query = """
[out:json];
(node"""+query+"""(around:"""+radious+""","""+lat+""","""+lon+""");
way"""+query+"""(around:"""+radious+""","""+lat+""","""+lon+""");
relation"""+query+"""(around:"""+radious+""","""+lat+""","""+lon+""");
);
out geom;
>;
out skel qt;
"""
response = requests.get(overpass_url,params={'data': overpass_query})
data = response.json()
data = osm2geojson.json2geojson(data)
json.dump(data,open(output,'w' ))

Below and example of the error that I get when I run the query "route" = "subway"
Failed to convert relation to shape {'ref': 8095058, 'role': 'platform', 'type': 'relation'}

Any thoughts?


texts: Openlayers 3 polygon not rendering on map

texts:
I've set up a small example of an OpenLayers 3 map, and I am trying to add a custom polygon to the map. If I use the data from the example, it works correctly, the polygon renders, and it is shown on the map, but when I try to add another custom polygon, it is not showing up on the map at all. Here's the code I use:

function initMap(map) {

 var IMap=new Map(map, [50.5973695, 6.200000]);
 var polygons=jQuery("#"+map).data('polygons');
 console.log(polygons);
 jQuery.each(polygons, function(index, value) {
     IMap.addPolygon(index, value);
 });

 IMap.addPolygon('test', [[-12.55,-21.12], [-44,-55], [-88,75]]);
}

Function Map works correctly, as the map is rendered, and centered to the specific position. Function addPolygon also works, because the polygon which I add manually with the IMap.addPolygon line renders on the map. Here is the Map::addPolygon function:

Map.prototype.addPolygon = function(name, coordinates) {

 var _coordinates=[], i;
 for (i=0; i<coordinates.length ;i++ ) {
     console.log(coordinates[i]);
     _coordinates.push(ol.proj.transform(coordinates[i], 'EPSG:4326', 'EPSG:3857'));
 }

 console.log("Creating polygon", _coordinates);
 var Polygon = new ol.geom.Polygon([_coordinates]);
 var feature = new ol.Feature({
     name: name,
     geometry: Polygon
 });
 this.polygons.push(feature);

 feature.setStyle(polyStyle);

 this.vectorSource.addFeature(feature);
}

Here is the two _coordinates which is logged out:

This one is for the dynamicaly creted polygon, which is not rendering:
Creating polygon Array [ Array[2], Array[2], Array[2], Array[2], Array[2], Array[2], Array[2], Array[2], Array[2], Array[2], 9 more… ]
Inner array example from above array: Array [ 50.5982266, 6.2467439 ]

This one is for the manually added polygon, which renders correclty:
Creating polygon Array [ Array[2], Array[2], Array[2] ]
Inner array example from above array: Array [ -13.12, -21.76 ]

What could be the problem?

NOTE: This question was also posted on stackoverflow.com.


texts: Set size of map in PageLayoutControl to print ScaleBar underneath, not within map

texts:
For printing we have a PageLayoutControl where the user can select the map extent, he wants to print. I like the map to fill only a square in the middle of the page so that the title, the ScaleBar and some additional texts appear on top and underneath the map. By now the map fills the hole page. How can I resize the map?

I tried setting the envelopes x and y values, but with no success:

public void SizeMap(IPageLayout pageLayout, IActiveView activeView)
{
  IMap map = activeView.FocusMap;
  IGraphicsContainer graphicsContainer = pageLayout as IGraphicsContainer;
  IFrameElement frameElement = graphicsContainer.FindFrame(map);
  IMapFrame mapFrame = (IMapFrame)frameElement;

  IElement mainMap_Element = mapFrame as IElement;
  IGeometry geometry = mainMap_Element.Geometry;
  IEnvelope mainMap_Envelope = geometry.Envelope;

  mainMap_Envelope.YMin = 5; // only example values
  mainMap_Envelope.YMax = 12;

  activeView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);
}

Is this the right approach?


texts: Exporting slope raster for entire country in GEE

texts:
I have loaded SRTM (30m resolution) DEM data through Google Earth Engine and have determined the slope for the entire country of Peru.
When I export the image out into Google Drive, I get the following error message- "Export region contains no valid (un-masked) pixels".

How can I export out slope rasters (tifs mostly) via GEE for entire countries?
My code is available here.

https://code.earthengine.google.com/2f490881ca64a744041c0cc702a24a95

I have also changed the max number of pixels that can be exported, but I still hit the same problem.


texts: ArcMap reports that PostGIS polygon data is self-interescting

texts:
I have an ESRI File Geodatabase of data that I upload to PostGIS using FME. Seeing as I have to have the data in PostGIS for other purposes, it would be great to view the data in ArcMap straight from the PostGIS database, rather than having to use the file geodatabase. I dont need to edit the data in ArcMap.

If I open the data from the ESRI File Geodatabase in ArcMap or QGIS, the data displays fine.

If I load the PostGIS version in QGIS, it displays fine.

However, when I open the PostGIS version in ArcMap, I get the following error:

  One or more layers failed to draw: Linestring or poly boundary is self-intersecting

In my FME workspace I am using a GeometryValidator transformer to check that the geometries are correct - I use the tests Fails OGC Simple and Fails OGC Valid and write only the valid geometries to PostGIS. At the moment I am only loading the valid geometries to PostGIS (I will deal with the invalid ones once I get get ArcMap to display the valid ones).

I have double checked that the data is valid in PostGIS by checking the data with ST_IsValid().

I am using ArcGIS 10.2.1 and what I believe is a supported version of PostgreSQL/PostGIS (PostgreSQL 9.2.4 64-bit and PostGIS 2.0.3). I also have the dll files on the client machine and have no problems connecting the the database.

I have also tried with a newer version of both PostgreSQL and PostGIS (which are not supported by ESRI) and they dont work either.

Any ideas?


texts: Getting India map on QGIS canvas using PyQGIS?

texts:
Can I get India map extent in one click from OSM tile on canvas using PyQgis?

This is my code snippet.

canvas = iface.mapCanvas()
urlWithParams = 'type=xyz&url=https://a.tile.openstreetmap.org/{z}/{x}/{y}.png&crs=EPSG3857'
rlayer = QgsRasterLayer(urlWithParams, 'MapMyIndia', 'wms')
QgsProject.instance().addMapLayer(rlayer)
rect = QgsRectangle(7589389.42, 889589.56, 10842803.55, 4231219.37)
canvas.setExtent(rect)
canvas.update()


texts: Open Data Cube: Unable to connect to database within opendatacube Conda environment

texts:
I followed the steps listed at https://cubes.readthedocs.io/en/v1.0.1/install.html
and did the following:

Installed Python

Installed PostgreSQL
ssudo sh -c 'echo "deb http://apt.postgresql.org/pub/repos/apt $(lsb_release -cs)-pgdg main" > /etc/apt/sources.list.d/pgdg.list'
wget --quiet -O - https://www.postgresql.org/media/keys/ACCC4CF8.asc | sudo apt-key add -
sudo apt-get -y install postgresql

Installed Miniconda and created a conda datacube environment
conda config --add channels conda-forge
conda create --name cubeenv python=3.6 datacube
source activate cubeenv

datacube system check

Gives me an error stating

psql: error: could not connect to server: could not connect to server: No such file or directory
Is the server running locally and accepting
connections on Unix domain socket "/tmp/.s.PGSQL.5432"?

Problem is that PostgreSQL is accessible outside the environment but not getting connected once the environment is activated.
A screenshot of the same is attached.


texts:
Still pretty new to GIS, but Ive been working on assembling a map with ~20 vector layers over the past couple of weeks. Its finally to a point where Im ready to export. Ive been experimenting with various export methods for a few days now, and I cant seem to get what Im looking for.

Map info -
its massive, covers over a million acres.
my goal is to have it printed on a banner, 10 feet wide, by however tall makes sense to preserve the aspect ratio
For prints of this size, DPI needs to be high, ~300 or more, from what Ive read
Scale - Im trying to keep it around 1:24000, or even 1:20000

I have been trying for either a PDF or a TIF file.

When I use the exact settings I want, the export fails with one of those errors about it being too big, cant allocate enough memory, ect.

This happens for both the Project -> export, and the print composer methods.

When I tweak the settings, compromising on either scale, or dpi, I eventually get a successful output, kind of. The files it creates look clean, open up fine, but they have part of map cut off, always the far right edge. And sometimes there is even a strange black block that appears in the area that was cutoff. It seems as though the exporting process works left to right across my map, and fails at some point, but doesnt report it as a failure.

Below is an example of what Im talking about, Ive stacked 3 different export attempts, with my full QGIS map on the top (please ignore the label mess):

Any suggestions for what else I might try?


texts: How to configure grouped layers in the LayerList widget with feature layers

texts:
This applies to the ArcGIS Server JavaScript API v4.x.
This LayerList widget example code contains grouped layers, whereby the San Francisco parent checkbox also controls its child layers:

The example code uses a WebScene document as the source for the layers, but in my case I'm manually creating the layers as FeatureLayer objects, and adding them to a Map object, using the pattern:
const myMap = new Map({basemap: "streets-vector"});
const layer1 = new FeatureLayer({url: "https://sampleserver6.arcgisonline.com/arcgis/rest/services/USA/MapServer/0"});
const layer2 = new FeatureLayer({url: "https://sampleserver6.arcgisonline.com/arcgis/rest/services/USA/MapServer/1"});
myMap.addMany([layer1,layer2]);

Using this pattern, is it possible to create layer groups/hierarchies, ie create a parent checkbox for layer1 and layer2?


texts: Application Extension vs. Add-in Application Extension

texts:
Background:
I am tasked with rewriting a 9.3 ArcMap solution based upon a VB6 implementation using an application extension and several forms.  The interface to ArcMap is that this solution needs to find a running instance of ArcMap with a specific map in it and then selects features in a layer in that map.  This will be redone for ArcMap 10.2 using C# and VS2012.  I am new to ArGIS stuff and in someways new to .NET/COM development.  Yes, I know I am in deep. :)

My Question:
I have read the documentation till my head is spinning - such a mixture of versions and FWs between both ESRI and MS.  I cannot decide upon the proper approach.  I feel that I should use a COM based application extension as the heart of the re-write, but I have also seen mentioned an Add-in application extension.  What I don't know is, can I use an Add-in application extension instead of a full blown COM based one?  Do they have the same abilities, or as in Add-in application extension more limited than a COM one? I have seen information that points to add-ins not being as capable as COM based solutions but this seems to be lessening with each release.  To be honest, I don't even know if an Add-in application extension is the same kind of thing as a "traditional" application extension - and I can't find that out from the docks, which frankly, suck.

Additionally, if I AM able to go with an Add-in application extension, how capable are the dockable forms that apparently are part of Add-in land?  In the VB6 Forms that I am replacing, database work using ADO to Oracle is being done.  Would I still be able to do that with these dockable windows?


texts: I would like to know the use cases of Esri CityEngine, and what distinguish this tool from the other standard GIS like ArcMap and QGIS? and in which fields could be used ?


texts: Openlayers – Empty spaces when resizing map in height (opera and safari)

texts:
The following test site is supposed to let the map extend to the maximum available space and to ensure a minimal height of 300px.
This works perfectly in firefox.
But in opera and safari on all apple devices (Mac, IPhone, IPad), after resizing the window a few times will create empty areas (screenshot).

Js fiddle

Jemand musste Josef K. verleumdet haben, denn ohne dass er etwas Böses getan hätte, wurde er eines Morgens verhaftet. »Wie ein Hund!« sagte er, es war, als sollte die Scham ihn überleben. Als Gregor Samsa eines Morgens aus unruhigen Träumen erwachte, fand er sich in seinem Bett zu einem ungeheueren Ungeziefer verwandelt. Und es war ihnen wie eine Bestätigung ihrer neuen Träume und guten Absichten, als am Ziele ihrer Fahrt die Tochter als erste sich erhob und ihren jungen Körper dehnte. »Es ist ein eigentümlicher Apparat«, sagte der Offizier zu dem Forschungsreisenden und überblickte mit einem gewissermaßen bewundernden Blick den

Ps:
At the time of writing the OSM map server seems to be down.
(Working again)


texts: Changing QGIS map from ' bird's eye / aerial view' to normal view?

texts:
While working in a map, the map changed in to bird's eye view'. I've managed to get the base layer back into 'normal' view by changing the coordinate reference. The other layers are still in aerial view and on a different 'location' even when I change their cooordinate reference.

Any idea how I can get it all layers back to normal view and on top of the base layer?


texts: Polygon Volume (ArcMap 10.7) calculates a volume of zero?

texts:
I'm calculating partial TIN volumes using Polygon Volume in ArcMap 10.7. I've done this hundreds of times using arcpy and for the most part, it works fine (Height Field = MAX; Reference Plane = BELOW). However, sometimes the routine calculates a volume and a surface area of zero, respectively, without any error messages.
My data can be downloaded here. I have provided one example of a TIN and a shapefile that works, and one example of a combination that doesn't. I cannot figure out the difference between the two.
EDIT:
Following @Vince, I've used Check Geometry on the shapefile (as it does not work with TINs), which gave me the error notification:
Could not find spatial index.
When Repair Geometry did nothing to fix this, I simply ran Add spatial index.
After that, the Check Geometry routine gave me no further errors - however, I still get an area and volume of zero when I use Polygon Volume.
When I use other shapefiles (that overlap with the TIN to a larger or smaller extent) the volume/area calculation runs fine. So, I suspect the error has to lie within the shapefile - but I have no idea what caused it.


texts: How to configure Geowebcache as standalone application?

texts:
I have a cluster of Geoserver and I am trying to use Geowebcache as a standalone application.

I downloaded the war file from this website:
https://sourceforge.net/projects/geowebcache/files/geowebcache/1.12.2/

I unzipped it inside a separate tomcat, restarted tomcat, modified the geowebcache.xml so it links to Geoserver layers and it worked fine (I can see my layers in Geowebcache demo page).

Now, how can I use this layers directly from Geoserver? For example, when a client makes a request, Geoserver will verify if the layer exists in the cache before creating it.


texts: What Information Does the Blue Band Provide?

texts:
I am playing around with ERDAS and was implementing image difference on two rasters separated by interval of time. I was using a blue band raster and I created a default threshold for percentage change to be 20% and also did a standard deviation of 2.0.

The image gives two completely different outputs for the two types of measuring the difference, using the percentage change and standard deviation. This is the image.

The green describes the lower extreme, describing the reduction in change and the red describes the higher extreme, indicating the increase in change. Though I am not able to understand the difference between the standard deviation output and the percentage output and what is the significance of using one over the other?

Hope I was able to explain the question well. In the image the upper one is the percentage change and the lower one is the standard deviation.


texts: PostGIS Rasters not aligned even after calling ST_Resample() or ST_SnapToGrid()

texts:
I am trying to compute certain statistics from an elevation (DEM) raster column but at certain points the raster tiles are not aligned, and continue to not be aligned even though ST_SnapToGrid() or ST_Resample() are called.

with
    buffer (buf) as
        (values (ST_Transform(ST_Buffer(
          ST_SetSRID(ST_Point(-66.0005555554938, 18.0005555553954),
               4269)::geography, 30)::geometry, 4269))),
    raster as
        (select
             ST_Union(ST_SnapToGrid(rast,
                -66.0005555554938, 18.0005555553954))
             as rast
         from raster_table, buffer where ST_Intersects(rast, buf))
select
    rast From raster;

with a similar approach using ST_Resample and another CTE to use the first raster returned from the intersection as a sort of reference. I have tried different resample methods and different values for maxerr, include 0, 1.0, and 100, but none of these prevent the error SQL Error [XX000]: ERROR: rt_raster_from_two_rasters: The two rasters provided do not have the same alignment from being raised. Is there something I'm missing or a better way to do this?


texts:
Has anyone had any success deploying the Ordnance Survey's free datasets on its open sourced OSGBWebMapTools framework using GeoServer 2.0.2?

I've been able to set up the datasets on GeoServer 1.7.6 and would prefer for it to use the out of the box caching that GeoServer 2.x supports.


texts: Using QGIS Atlas outside of Print Composer?

texts:
I wonder if it is already possible to use the Atlas Panel outside of Print Composer like it would be in ArcGIS with Data Driven Pages. I often use Data Driven Pages to screen through an entire area in order to edit/fix the results of a vectorized classification. Basically I use an index table to check every box of my classification. However, since I'd like to migrate to a free software, I would love to be able to do the same on QGIS.

I wonder if this is already possible of if there is any plugin that would make it possible to use an index table as an auxiliary file to do the screening of an entire raster.


texts: Internet mapping UX design strategy: Separate desktop and mobile sites, or one site fits all?

texts:
I searched GIS.SE for "Desktop vs. Mobile" and didn't get anything helpful, so I guess I'll be the one to ask.

Question: With respect to Internet mapping applications, is it better to have separate sites for desktop and mobile, or to pursue one highly adaptable site (Responsive Web Design) that attempts a best-fit across all devices?

Background: A quick review of this problem reveals this article (Nielsen, April 10 2012) where the author advocates separate sites for desktop and mobile, as well as this rather strong rebuttal (Lawson, April 19 2012) encouraging development of one well-planned site capable of serving a variety of devices. Other articles discuss various pros and cons, but ultimately remain neutral. [1] [2] [3]

My Own Experience: My first attempt at this was the one-site-fits-all approach, and while it renders properly on most mobile devices, I find it lacking with respect to the total experience. I want better use of screen area. I want bigger control surfaces to prevent "fat fingering" consequences. And frankly, unlike Lawson, I think our mobile users probably do have different goals than our desktop users. I doubt our mobile users care about our 2-foot elevation contours, but I bet they commonly want as much parcel data as they can get, preferably in a layout that's convenient on mobile.

At the same time, maintaining two separate web mapping clients seems like a bad idea. ..then again, maybe that's a false economy? Perhaps two simpler, more focused sites would be easier to maintain than one burgeoning, more complex site?

Denouement: The Lawson article suggests there will always be edge cases favoring multiple sites over responsive design, and I'm starting to believe Internet mapping is one of the edge cases. So I'm asking the community to provide insight. Which of these strategies do you prefer and why? I figure the better strategy should organically up-vote its way to the top, providing a kind of de facto authentication of the approach that's better suited for Internet mapping applications.


texts: Changing coordinate systems in ArcGIS Desktop?

texts:
I am working on a map in ArcGIS ArcMap v10.2 with coordinate system GDA94 (UTM). I have just added a layer with points in latitude and longitude (WGS84). I specified that the layer was WGS84 and the points appear to be in the correct position on the map.

Do I need to use ArcToolbox (Data Management - Projections and Transformations - Feature - Project) to change this coordinate system for the final map, or is it okay as is?


texts: refreshLayerSymbology raster layer legend

texts:
When I modify raster colors with PyQGIS, the legend in the layout layer is not update (color and value). The function refreshLayerSymbology() in my code don't do the job.

I had try this :

# BBOX
extent = box
# acces données
provider = self.MNT.dataProvider()
stats = provider.bandStatistics(1, QgsRasterBandStats.All, extent, 0)

if (stats.minimumValue < 0):
    min = 0

else:
    min = stats.minimumValue
max = stats.maximumValue
range1 = max - min
add = range1 // 2
interval = min + add

colDic = {'red': '#ff0000', 'yellow': '#ffff00', 'blue': '#0000ff'}

valueList = [min, interval, max]

lst = [QgsColorRampShader.ColorRampItem(valueList[0], QColor(colDic['blue'])),
       QgsColorRampShader.ColorRampItem(valueList[1], QColor(colDic['yellow'])),
       QgsColorRampShader.ColorRampItem(valueList[2], QColor(colDic['red']))]

myRasterShader = QgsRasterShader()
myColorRamp = QgsColorRampShader()

myColorRamp.setColorRampItemList(lst)
myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
myRasterShader.setRasterShaderFunction(myColorRamp)

myPseudoRenderer = QgsSingleBandPseudoColorRenderer(self.MNT.dataProvider(),
                                                    self.MNT.type(),
                                                    myRasterShader)
# Add
self.MNT.setRenderer(myPseudoRenderer)

# Refresh the canvas and the legend
self.MNT.triggerRepaint()
#self.iface.mapCanvas().refresh()
self.iface.legendInterface().refreshLayerSymbology(self.MNT)

However when I use the function refreshLayerSymbology() with grayscale color, the function work.

The code I used is :

myGrayRenderer = QgsSingleBandGrayRenderer(self.MNT.dataProvider(), 1)

self.MNT.setRenderer(myGrayRenderer)

renderer = self.MNT.renderer()
# acces données
provider = self.MNT.dataProvider()

layer_extent = self.MNT.extent()
uses_band = renderer.usesBands()

myType = renderer.dataType(uses_band[0])

stats = provider.bandStatistics(uses_band[0],
                                QgsRasterBandStats.All,
                                layer_extent,
                                0)

myEnhancement = QgsContrastEnhancement(myType)
contrast_enhancement = QgsContrastEnhancement.StretchToMinimumMaximum

myEnhancement.setContrastEnhancementAlgorithm(contrast_enhancement, True)
myEnhancement.setMinimumValue(stats.minimumValue)
myEnhancement.setMaximumValue(stats.maximumValue)

self.MNT.renderer().setContrastEnhancement(myEnhancement)

# Refresh the canvas and the legend
self.MNT.triggerRepaint()
self.iface.legendInterface().refreshLayerSymbology(self.MNT)


texts: Cannot open Shapefile with any Python libraries, but it will open in QGIS/ArcGIS

texts:
I have a shapefile that I cannot open in any open source Python libraries, including GDAL, Pyshp, Fiona, and Geopandas.  The error message: Invalid length for entity 9164 is a GDAL error message and can be found in the source code: https://fossies.org/linux/gdal/ogr/ogrsf_frmts/shape/shpopen.c at lines 871 and 2121.  Both sections indicate that the error occurs during the reading/handling of offsets in the .shx file.  Here is the comment to the first section that the error message occurs at: Read the .shx file to get the offsets to each record in the .shp file

The really frustrating aspect of this issue is that the file can be viewed and edited in QGIS/ArcGIS with no issues.

If I open the shapefile in QGIS and then save it with a different name, I can open the new shapefile with all python libraries.  Upon further investigation, if I use the new .shx file with the original set of shp files, the file also opens in Python.  So somehow this is a corrupted .shx file.  Any info would be helpful as I have my processes automated.


texts: What is the default unit of Zonal Geometry Table and how to change it?

texts:
I want to summarize the area of each region group of my result. I sort them out based on the value (score that I granted them). Now I need to pick the best 15% land from the targeted area based on the score.

My method now is to list the values of all the region groups and I add up the area of first 10-15 group maybe until it fills the goal. However the info about the land I have is in acre, which is numerically smaller and easier to calculate. And the result of Zonal Geometry Table is largely and clearly not sharing the same unit as shown below

What is the unit of the AREA column and is there anyway to change that?


texts: Creating Store in GeoServer workspace using REST API?

texts:
I'm building a PHP class to interact with my GeoServer. I'm following the limited Documentation, this GeoServer REST API with PHP post and also this great advice. I've successfully created/deleted a Workspace.

But when I try to upload a Shapefile it fails because no dataStore endpoint exists in my Workspace (I get a 404 response from cURL). So, I guess I should create a Store before trying to upload a Shapefile.

I've failed at finding any information regarding this issue (both here at GIS SE or Googling) but I've successfully created a Store with some RevEng:

I've checked the JSON for my GeoServer's default Store:

{
  "dataStore": {
    "name": "cntry_shp",
    "type": "Directory of spatial files (shapefiles)",
    "enabled": true,
    "workspace": {
      "name": "opengeo",
      "href": "http:\/\/0.0.0.0:8080\/geoserver\/rest\/workspaces\/opengeo.json"
    },
    "connectionParameters": {
      "entry": [
        {
          "@key": "memory mapped buffer",
          "$": "false"
        },
        {
          "@key": "timezone",
          "$": "America\/Vancouver"
        },
        {
          "@key": "fstype",
          "$": "shape-ng"
        },
        {
          "@key": "create spatial index",
          "$": "true"
        },
        {
          "@key": "charset",
          "$": "ISO-8859-1"
        },
        {
          "@key": "filetype",
          "$": "shapefile"
        },
        {
          "@key": "cache and reuse memory maps",
          "$": "true"
        },
        {
          "@key": "enable spatial index",
          "$": "true"
        },
        {
          "@key": "url",
          "$": "file:workspaces\/opengeo\/data"
        },
        {
          "@key": "namespace",
          "$": "http:\/\/opengeo.org"
        }
      ]
    },
    "_default": false,
    "featureTypes": "http:\/\/0.0.0.0:8080\/geoserver\/rest\/workspaces\/opengeo\/datastores\/cntry_shp\/featuretypes.json"
  }
}

I copied it, replacing its name and workspace, and inserted it in my system with a cURL call at the endpoint:

 rest/workspaces/MY_WORKSPACE/datastores

Then I tried to upload a Shapefile to it, and it failed. Anyway, I don't think this is the right way to do this.

Do you know which is the right way or, at least, do you know some examples or documentation?


texts: Batch unsupervised classification with ERDAS 2011

texts:
I work with Erdas Imagine 2011. I would like to do an unsupervised classification and vectorization on 100 Landsat images (.tif).

CC_2014_1.tif
CC_2014_2.tif
CC_2014_3.tif
etc

However, I am new in Erdas. Please can you indicate me few materials or scripts which can help me to do the unsupervised classification and the vectorization on my images.


texts:
I have ArcPad 10 installed on a Trimble Juno.  I tried to use a geodatabase with a feature class set up with domains and subtypes for field data collection, but, I learned quickly that ArcPad10 doesn't recognize a geodatabase... only shapefiles.

I understand that an *.apl file (part of a shapefile package) is used in arcpad10 to constrain data editing, or create forms.

Do you need ArcStudio to create an apl file?

What is the structure of an apl file?

Can you write your own?

I am running ArcGIS Desktop 9.3, and I see the ArcPad extension but I am assuming you need ArcPad installed on your desktop. I only have it installed on my Juno Trimble. I will not have access to the ArcPad software right away.


texts: Making routable network topology using pgRouting?

texts:
We are trying to make a network topology for pgRouting from our own data. But there is a problem we could not resolve; our road data look like this;

The red edge is a single geometry, so we can not make a route between these nodes.
How can we convert this data to be able to make routing work?


texts: On-the-fly reprojection to a "world" projection causes lat/lon raster to become invisible

texts:
I have a global elevation raster in latitude and longitude. I want to use it as the basemap in a world projection system (e.g., Mollweide, Robinson).

When I set the coordinate system to Mollweide (or Robinson) and enable on-the-fly projection, the global topography raster map disappears until I zoom in far enough that there is no "no data" area left in the display. In other words, I need to make sure that the display doesn't include any of the curved map projection edges.

I have tried using gdalwarp outside of QGIS as well, and then importing the already-projected raster, using gdalwarp:

gdalwarp -s_SRS 'EPSG:4326' -t_srs '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs' coarse_6am_topo_021000.tif out.tif

However, this projected map cuts off the north polar region (not sure why this is happening either), so I am not sure what to do -- it would be nice if the on-the-fly reprojection worked!


texts:
Environment: Windows 7, ESRI ArcGIS Pro 2.1, R 3.5.1 with reticulate v1.10.

ArcGIS Pro comes with a pre-configured python 3 conda environment that I'm trying to access in R. Related questions here and at stackoverflow all seem to focus on a similar situation with ArcMap 10.x, but that uses python 2.7 and doesn't use conda, so they haven't been much help.

If I'm reading the reticulate documentation correctly, the following should work:

library(reticulate)

arcpy3_dir <- 'C:/Program Files/ArcGIS/Pro/bin/Python'

use_condaenv(condaenv = file.path(arcpy3_dir, 'envs', 'arcgispro-py3'),
             conda    = file.path(arcpy3_dir, 'Scripts' ,'conda.exe'),
             required = TRUE)

However, I'm getting an 'unable to locate conda environment' error.

There are some .bat files in the same directory as conda.exe that seem to do a lot of temporary messing about with system paths (more commentary here for an older version of ArcGIS Pro), but I don't know enough about shell scripts to clearly understand what they're doing. I think I need to set some environment variables with Sys.setenv(), but I'm not sure which. I tried

Sys.setenv('CONDA_NEW_ENV' = 'C:/Program Files/ArcGIS/Pro/bin/Python/envs/arcgispro-py3',
           'CONDA_EXE'     = file.path(arcpy3_dir, 'Scripts' ,'conda.exe'),
           'CONDA_SCRIPTS' = file.path(arcpy3_dir, 'Scripts'))

but that doesn't appear to be sufficient.

Has anyone managed to successfully access the ArcPro conda environment in R, using reticulate?


texts: Loading mod_spatialite using sqlite3 for Python in Windows

texts:
I have been successfully using the mod_spatialite extension for sqlite in Linux for a while now, but a client wants an application that relies on mod_spatialite to work in Windows.

I have extracted the required dlls to a folder that I have added to the Windows PATH variable. This lets me start a standalone version of sqlite3.exe and load mod_spatialite with SELECT load_extension('mod_spatialite');. However, when trying to load this extension in Python using the sqlite3 module, I get an error:

  sqlite3.OperationalError: The specified procedure could not be found

In linux, I can simply use SELECT load_extension('/path/to/mod_spatialite');, but sqlite3 gives the same 'specified procedure could not be found' error.

Am I missing something here?


texts: How to change the view from layout view to data view using c# code in ArcMap?

texts: How to change the view from layout view to data view using c# code?


texts: Using arcpy.SetParameter() with Python Script Tool published as Geoprocessing Service?

texts:
I have a script that needs to be published to ArcGIS Server as a geoprocessing service.  This script computes a set of records in table form.  I need to set an output parameter so that it returns the table (recordset).

So, to make sure I have a table (not a feature class), I am using arcpy.TableToTable_conversion(), outputting the table to the in_memory workspace.  So, I know I have a table.  I know that it has records, as I can perform an arcpy.GetCount and get the number of records.

I have the output parameter defined.  Here's what it looks like:

arcpy.SetParameter(2, theTable)

No matter what I try this does not work.  In my toolbox, the output parameter is set as a RecordSet.  I've tried supplying a recordset object, the path to the in-memory feature class, and many other ideas.

How do you properly handle outputting a table to a recordset output using SetParameter for ArcGIS server?


texts:
I need to know if in QGIS 3.16 there is a method to stretch a label all along line length.
My layer has lines with different length and I want to know if is possible to stretch automatically labels all along lines length, preferably using letter spacing. I've never used rule-based label and I don't know how it works, maybe I should use it to do what I want?


texts: QGIS raster merge results in a black rectangle

texts:
I have 4 SRTM raster like this:

I want to merge them into 1, but when I do this, something like this happens (I only tried it with two of them here):

I also tried changing the min and max values that are here both 0 to the min and max values from the two used rasters but it doesn't change.
Do you have any idea what's the problem here and what I can do?

Here's the Log file:

  QGIS version: 3.8.0-Zanzibar QGIS code revision: 11aff65f10 Qt
  version: 5.11.2 GDAL version: 2.4.1 GEOS version: 3.7.2-CAPI-1.11.0
  b55d2125 PROJ version: Rel. 5.2.0, September 15th, 2018 Processing
  algorithm… Algorithm 'Merge' starting… Input parameters: { 'DATA_TYPE'
  : 5, 'INPUT' : ['C:/Users/.../srtm/eu_dem_v11_E30N20/eu_dem_v11_E30N20.TIF','C:/Users/.../srtm/eu_dem_v11_E30N30/eu_dem_v11_E30N30.TIF'],
  'NODATA_INPUT' : None, 'NODATA_OUTPUT' : None, 'OPTIONS' : '',
  'OUTPUT' : 'TEMPORARY_OUTPUT', 'PCT' : False, 'SEPARATE' : False }

  GDAL command: python3 -m gdal_merge -ot Float32 -of GTiff -o
  C:/Users/.../AppData/Local/Temp/processing_0bf9ee3ec2f548829744a6bd244c96a2/60f16fc4e2394ab496e96b39dd323a73/OUTPUT.tif
  --optfile C:/Users/.../AppData/Local/Temp/processing_0bf9ee3ec2f548829744a6bd244c96a2/1eb95beda25545199dda413988a6c112/mergeInputFiles.txt
  GDAL command output: 0Traceback (most recent call last):

  File "C:\PROGRA~1\QGIS3~1.8\apps\Python37\lib\runpy.py", line 193, in
  _run_module_as_main

  "main", mod_spec)

  File "C:\PROGRA~1\QGIS3~1.8\apps\Python37\lib\runpy.py", line 85, in
  _run_code

  exec(code, run_globals)

  File "C:\PROGRA~1\QGIS3~1.8\bin\gdal_merge.py", line 611, in

  sys.exit(main())

  File "C:\PROGRA~1\QGIS3~1.8\bin\gdal_merge.py", line 596, in main

  fi.copy_into(t_fh, band, band, nodata)

  File "C:\PROGRA~1\QGIS3~1.8\bin\gdal_merge.py", line 343, in copy_into

  nodata_arg)

  File "C:\PROGRA~1\QGIS3~1.8\bin\gdal_merge.py", line 124, in
  raster_copy

  m_band)

  File "C:\PROGRA~1\QGIS3~1.8\bin\gdal_merge.py", line 186, in
  raster_copy_with_mask

  to_write = Numeric.choose(mask_test, (data_src, data_dst))

  File
  "C:\PROGRA~1\QGIS3~1.8\apps\Python37\lib\site-packages\numpy\core\fromnumeric.py",
  line 422, in choose

  return _wrapfunc(a, 'choose', choices, out=out, mode=mode)

  File
  "C:\PROGRA~1\QGIS3~1.8\apps\Python37\lib\site-packages\numpy\core\fromnumeric.py",
  line 56, in _wrapfunc

  return getattr(obj, method)(*args, **kwds)

  MemoryError

  Execution completed in 941.61 seconds


texts: Null values in Geographically Weighted Regression Analysis

texts:
Some of my explanatory variables have a few null values for certain polygons.

Can Geographically Weighted Regression Analysis handle these, or should I remove the whole polygon and all data (dependent and independent variables for that particular polygon)?


texts: Is there a way to make a leaflet map's legend draggable? Like a library, or something?


texts: How to add attribute table to new shapefile using OGR?

texts:
i have made a new shapefile from the original shapefile , and now i want to copy the attribute table from the original shapefile to the new select shapefile , but i don't know how to do that?

import ogr,os,sys,ogr

Datasource_Parks = 'ArcCR500_v32/ChranenaUzemi.shp'

Direction_Result_Shp = 'ArcCR500_v32/Result.shp'
Direction_Result_Dbf="ArcCR500_v32/Table.dbf"

Open_Parks = ogr.Open(Datasource_Parks)
layer_park = Open_Parks.GetLayer(0)

SQL = "NAZEV_ASCI = 'Krkonossky narodni park'"

layer_park.SetAttributeFilter(SQL)

Number_Of_Park = layer_park.GetFeatureCount()

Driver = ogr.GetDriverByName("ESRI shapefile")

if os.path.exists(Direction_Result_Shp):
    Driver.DeleteDataSource(Direction_Result_Shp)

Dataset = Driver.CreateDataSource(Direction_Result_Shp)

New_Layer = Dataset.CreateLayer(Direction_Result_Shp, geom_type= ogr.wkbPolygon)

number = 0
for i in range(0, Number_Of_Park):
    H = layer_park.GetNextFeature()
    H_geometry = H.GetGeometryRef()

    New_Layer.CreateFeature(H)

Dataset.Destroy()


texts: Python Egg Extraction Error When Trying to Use pyproj in Pyspark

texts:
I have pyproj installed in my conda environment on my local machine (MacBook Pro).  Everything works properly.  I am able to use pyproj within the pyspark installation on my local machine.  However, I mainly work on a Linux machine I SSH to and I have shipped my conda environment from that Linux machine to an HDFS cluster so that I can use pyproj (and pyspark) on a large dataset.

When I try to run the same rdd.map() code on that HDFS cluster in which some pyproj stuff is inside the .map(), I get the following error:

ExtractionError: (ExtractionError("Can't extract file(s) to egg cache

The following error occurred while trying to extract file(s) to the Python egg
cache:

  [Errno 13] Permission denied: '/home/.cache'

The Python egg cache directory is currently set to:

  /home/.cache/Python-Eggs

Perhaps your account does not have write access to this directory?  You can
change the cache directory by setting the PYTHON_EGG_CACHE environment
variable to point to an accessible directory.
",).

The '/home/.cache/Python-Eggs' is a directory on the Linux machine I SSH to, not in HDFS, nor my local machine (MacBook Pro).

How can I debug this?


texts:
I am currently experimenting with the GeoServer ReST-API andy trying to update the configuration parameters of an ImageMosaic Coveragestore. I've been using Python 2.7 with the requests package as well as cURL, but after sending my POST-Request the layer isn't marked as active anymore in the GeoServer web interface.

I've saw this Message on the geoserver user mailing list and tried to reproduce this but without having success. No matter if I'm using Python or cURL, the POST-request does not update the specific entry in the parameters but overwrites all entries with only the one I specified in my request.

So here is my Python Code:

import requests

url = "http://MyGeoServer:8080/geoserver/rest/workspaces/Imaging_DOP/coveragestores/Geobasis_1410AB/coverages/Geobasis_1410AB.xml"
headers = {'Content-Type': 'text/xml'}
auth = ('user', 'password')

data = """
<coverage>
<parameters>
    <entry>
        <string>AllowMultithreading</string>
        <string>true</string>
    </entry>
    <enable>true</enable>
</parameters>
</coverage>
"""
r = requests.put(url, headers=headers, auth=auth, data=data)
print r.status_code

And here is the same request with cURL:

curl -u user:password -v -XPUT -H 'Content-type: text/xml' -d '<coverage><parameters><entry><string>AllowMultithreading</string><string>true</string></entry><enable>true</enable></parameters></coverage>' http://MyGeoServer:8080/geoserver/rest/workspaces/Imaging_DOP/coveragestores/Geobasis_1410AB/coverages/Geobasis_1410AB.xml

Does anybody know how to update specific parameters in CoverageStores without replacing all the entries, except the one specified in the request and without making the layer unusable?


texts: Solution for slow display of postgis-tables ? (qgis 1.8 + postgis 2.x)

texts:
I have the problem that my postgres-tables are displayed very slow when I scroll them in QGIS. I have no problem with other tested GIS (e.g. UDig, GVSig, ArcGIS). The number of objects isn't large (5000-10000) and the rendering is no problem.

I've observed that every action like scrolling causes queries in the database, looks like there's no cache - Is that true??? 'Can't believe that.

The connection to the database server runs through a fast Internet (and it works with other GIS). Does anyone have an idea how to solve the problem? I would love to work with QGIS, but the problem makes it impossible so far.

I appreciate your advices, regards Bob


texts: Search GeoNetwork metadata based on geonet:info:createdate or changedate

texts:
I want to retrieve metadata from the catalog with a request base on either changedate or createdate from geonet:info the additional information container inside the metadata
The following request returns the metadata base on the timestamp not the same as the one inside geonet:info.
<?xml version="1.0"?>
<csw:GetRecords
    xmlns:csw="http://www.opengis.net/cat/csw/2.0.2"
    service="CSW"
    version="2.0.2"
    resultType="results"
    outputSchema="csw:IsoRecord"
    maxRecords="%d"
    startPosition="%d">
    <csw:Query
        typeNames="gmd:MD_Metadata">
        <csw:ElementSetName>full</csw:ElementSetName>
        <csw:Constraint
            version="1.1.0">
            <Filter
                xmlns="http://www.opengis.net/ogc"
                xmlns:gml="http://www.opengis.net/gml">
                <PropertyIsGreaterThanOrEqualTo>
                    <PropertyName>createDate</PropertyName>
                    <Literal>%s</Literal>
                </PropertyIsGreaterThanOrEqualTo>
            </Filter>
        </csw:Constraint>
    </csw:Query>
</csw:GetRecords>


texts:
I'm looking for a method to change the initial view of the Qgis2threejs export. Is it very zoomed and I want less zoom when I open it.

Actually when I open the url is zommed:

But my idea is to have less initial zoom to see all the Italy.
I didn't find any option into the files to change it.

Any ideas?


texts: Remove all polygons except the main one from GeoJSON / Shapefile

texts:
I've got the coastline of Canada from OSM. For my work I only need the main landmass without any islands. Despite I was able to manually delete all the islands in JOSM / QGIS I need an automatic solution.

How to remove all polygons except the main one from a GeoJSON (or Shapefile) in commandline?

Even this simply doesn't work

ogr2ogr -f GeoJSON -dialect SQLITE -sql "SELECT * FROM Canada_1428125_AL2 order by ST_AREA(geometry) desc limit 1" out.geojson Canada_1428125_AL2.shp

because the multipolygon contains all islands. Is there any splitter function?


texts:
I have a raster loaded into PostGIS, but I do not know what bands it contains.  \d doesn't tell me.

How do I get a list of bands in the raster? (And anything else you feel I should know)


texts: QGIS crashes when adding vector layer from mysql or testing connection to mysql

texts:
For a while I have been adding layers from my mySql database and all went well.

Then I got crashes when adding new vector layers.

Then I only got crashes and it is not possible to add layers from the mySql database any more (crash) or create a connection, then test it (crash). Sporadically I can get to the list of the tables in the mySql database after creating a new connection but usually QGIS crashes.

I tried deinstalling and reinstalling QGIS. Did not help.

I am using the menu "layer" > "add vector layer" > "database" (in German, so maybe a little different).

The table I am trying to add is named "tblTeilPopMassnPolygon", contains an id "TPopMassnPolygonId" with Primary Key and Auto Increment. Plus a polygon field "TPopMassnPolygonPolygon" and a field named "TPopMassnGuid", varchar 40, to save GUID's in. There are 6000 datasets in the table. I had tested them with the polygon data before and it worked. So there is probably nothing wrong with the data itself.

I tried using shorter names for table and fields but that didn't help.

How can I get QGIS to connect to my mySql tables again?


texts: How to detect "delete" or "edit" events in popup menu created with Leaflet-draw-toolbar plugin?

texts:
In this example you may create a polygon and then click it and popup menu shows. There is three buttons: edit, delete and pick a color. How can I detect event and determine which of it - delete or edit - was executed in popup menu? I want to perform some action when one of this events is executed e.g. when a polygon is edited or deleted. The standart way for it:

map.on('draw:created', function (evt) {
        layer = evt.layer;
        // do something when polygon is created
});

map.on('draw:edited', function (evt) {
        layer = evt.layer;
        // do something when polygon is edited
});

map.on('draw:deleted', function (evt) {
        layer = evt.layer;
        // do something when polygon is deleted
});

But for events in popup menu it doesn't work.

So I edited file leaflet.draw-toolbar.js and added my custom event 'draw-popup:edited'

LeafletToolbar.EditAction.Popup.Edit = LeafletToolbar.ToolbarAction.extend({
...
    map.on('click', function () {
        shape.editing.disable();
        shape.fire('draw-popup:edited', shape);
    });
...
});

Now in my script file I'm able to detect 'draw-popup:edited' event

layer.on('draw-popup:edited', function() {
        var geom = this.toGeoJSON();
        console.log(geom);

});

But now it's fired multiple times when I click on a map and returns all shapes that were edited.
What is the right way to detect event and return the shape that was edited?


texts: Make a draw object option in Leaflet and that stores it in the webmap in a database or on a server

texts:
I would like to set up a simple web map which should be used as an information for those out in the field that needs to show something immediately. So I would just load the web map and get an overview.
The thing I would like to have is a map normal background map and then the user should be able to draw a polygon/point/line and that polygon should then be stored in a database or on the server.

For now I am going with Leaflet and QGIS.

Could someone tell me how it do this or provide links?


texts: Why do I get "Data source is invalid" error when calling method?

texts:
QGIS 2.14 on Windows 7. I'm trying to call a method in a python script from the QGIS python console. I've put the directory in which the scripts resides in the PYTHONPATH in both the windows system and user environment variable.

The python script (first_script.py) has a method that looks like this:

from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
def load_layer():
  wb = QgsVectorLayer('C:\Data\scripts\pyqbook\data\world_borders.shp', 'world_borders', 'ogr')
  QgsMapLayerRegistry.instance().addMapLayer(wb)

I call the script like this:

>>>import first_script
>>>first_script.load_layer()

and I get QGIS error: "Data source is invalid".
However, when I then manually enter the two lines of the method, it works fine and loads the layer into the canvas correctly:

>>>wb = QgsVectorLayer('C:\Data\scripts\pyqbook\data\world_borders.shp', 'world_borders', 'ogr')
>>>QgsMapLayerRegistry.instance().addMapLayer(wb)

Why can't I call the method from the QGIS python console directly without getting the "Data source is invalid" error? The data source is valid because I can call the individual lines of the method and get the correct result.


texts:
I'm using GEE in order to extract evapotranspiration (ET) data from MODIS (MOD16A2v6).

I want to extract only ET band in order to make some further process.
To do so, I apply some filters (like filterDate and a clip over my area of interest). Then I want also use the Quality Band information (bits) in order to create some mask function to keep only best quality data.

The problem is that when I want to load the image collection, after I applied my masks, GEE does not allow me to select only the ET band in a new variable.
I don't understand why. In my way of thinking, I would say that I load the entire Modis product collection with all its bands. Apply my filters and Masks on that collection (based on the QA band) and then select only the first product (band ET). But obviously it's not correct, I must be wrong somewhere.

Here is the part of my code related to that problem:

//*************** FUNCTIONS DEFINITIONS *******************

// ------------ Selection of data on the study area
var clip = function(img){
  return img.clip(subWet);
}; // détermine sur quelle zone on veut extraire les images modis. Ici, les 6 lagunes

// ------------- Filters

/* Create a function that select images on basis of QA information (stored in QA band).
Return only images according to specified bits.

  Function Arguments:
*   image - The QA Image to get bits from.
*   start - The first bit position, 0-based.
*   end   - The last bit position, inclusive.
*   name  - A name for the output image.
*/

var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
      pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// --- Mask out according to MODLAND_QC bits QA band (bit 0)
var mask1 = function(image) {
  // Select the QA band.
  var QA = image.select('ET_QC');
  // Get the internal_cloud_algorithm_flag bit.
  var internalCloud = getQABits(QA, 0, 1, 'Modland_QC');
  // Return an image masking out cloudy areas.
  return image.updateMask(internalCloud.eq(0));
};

// --- Mask out cloudy pixels - cloud state  from state QA band (bits 3 and 4)
var mask2= function(image) {
  // Select the QA band.
  var QA = image.select('ET_QC');
  // Get the internal_cloud_algorithm_flag bit.
  var cloud = getQABits(QA, 3, 4, 'cloud_state')
                        .expression("b(3) == 1 || b(3) == 2"); // sélection des deux bits que l'on dégage
  // Return an image masking out cloudy areas.
  return image.updateMask(cloud.not())//.internalCloud.eq(0));
};

// --- Mask out according to quality score -   from state QA band (bits 5to 7)
var mask3= function(image) {
  // Select the QA band.
  var QA = image.select('ET_QC');
  // Get the internal_cloud_algorithm_flag bit.
  var qscore = getQABits(QA, 5, 7, 'quality score')
                        .expression("b(5) == 2 || b(5) == 3 || b(5)==4"); // sélection des deux bits que l'on dégage
  // Return an image masking out cloudy areas.
  return image.updateMask(qscore.not())//.internalCloud.eq(0));
};

// ***************** CHARGEMENT DES IMAGES ****************

var dataset = ee.ImageCollection(MOD16A2v6)
var datafiltered=dataset
                  .map(clip) // définie pour notre zone clipée
                  .filterDate('2001-01-01', '2001-05-25')
                  .map(mask1)
                  .map(mask2)
                  .map(mask3)// function QA bits

var ET = datafiltered.select('ET') // sélection de la bande ET dans l'ensemble du produit modis


texts: Create diagonals in polygon with QGIS or ArcGIS

texts:
How is it possible to create the diagonals in/of polygons and to calculate the length of this diagonals afterwards.

It's for a lot of polygons so maybe a fast/automatic method would be nice.

I prefer the software QGIS or ArcGIS.

The polygons are created by following streets that cross each other.
So the shape is very weird, not a typically square.
Some polygons has 4 vertices, some only 3 vertices.

The green lines on the second figure shows the objective, it's to create this lines and calculate the length of this lines...

I am asking this question because of the following situation:
I digitalise this polygon on raster data of year 1950. Then i search the same ways in the raster data of year 2014. I try to digitalize the same polygon as in the year 1950.

But unfortunate the area (inside the ways) aren't perfect the same.
So i search something to compare the polygons of the same location of the two different years.
I was thinking to compare the area or perimeter, but isn't perfect.
so my next idea was to create the diagonals of the polygons and calculate the length.
If the difference in the length between two diagonals (of same polygon but different years) is to large, the two polygons aren't the same and so not possible to use in the further analysis...
maybe there is an easier way to compare the polygons?


texts:
I have an application where the user selects some "markers" from a google maps (all those coordinates' markers are stored into a postgresql database).

Right now, I can draw the route between all selected markers (using DirectionService) but when I have over 10 waypoints the route seems to be drawn in the user selected order (not in the "shortest route" order - I also used optimizeWaypoints parameter to true and split the route into subroutes of 8 points + start and stop points without success).
What I`m trying to do now: when the user selects a marker, I save his selection (marker coordinates) and then all selected coordinates are sent (via ajax)
 to the server.

The received coordinates are in a format something like this:

array((lat1,lng1), (lat2,lng2), ...(latn,lngn))

From here, I want to return back to the user all those coordinates ordered by "shortest path" and I don't know how to achieve this.

Any idea will be great.


texts: Pass temporary output between definitions / functions to use as input to next def

texts:
Lets say I have two functions:

class AnalyseDTM():

    def __init__(self):
        self.DTM = ('C:/Users/EH/QGis/Kvit_Tiny.tif')
        self.CreatePolygons()
        self.test()

    def CreatePolygons(self):

        translateFeatures = processing.runAndLoadResults("gdal:translate", {
        'INPUT':self.DTM,
        'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:5110'),
        'NODATA':0,
        'COPY_SUBDATASETS':False,
        'OPTIONS':'',
        'EXTRA':'',
        'DATA_TYPE':0,
        'OUTPUT': 'TEMPORARY_OUTPUT'})
        return

    !HERE!

        RemoveSub0values = processing.runAndLoadResults("qgis:rastercalculator", {
        'EXPRESSION':'((\"OUTPUT@1\">0)*\"OUTPUT@1\") / ((\"OUTPUT@1\">0)*1 + (\"OUTPUT@1\"<=0)*0)',
        'LAYERS': translateFeatures['OUTPUT'],
        'CELLSIZE':None,
        'EXTENT':None,
        'CRS':QgsCoordinateReferenceSystem('EPSG:32633'),
        'OUTPUT': 'TEMPORARY_OUTPUT' })

I've removed about 20 other processing functions as the above should hopefully explains what I'm looking for. So you see I have a def CreatePolygons, which currently has two processing. functions in it. This works fine, however I'm wondering if where I have written !HERE! I can separate the two functions. For example I could call it def test(self):

Therefore, I'm wondering if I can pass the 'OUTPUT' from translateFeatures and have it as the input for RemoveSub0values even when they are in separate def 's. I know that if I save it to file and have it in the __init__ then yes, I can, as this is part of a large model I don't want to save multiple outputs to file every time.

I have a fully working class, but its extremely slow, and I therefore wanted to test if grouping processing algorithms would speed this up, but I therefore need to be able to pass temporary outputs between functions.


texts: Is there any online or desktop based tool that can convert shapefile or feature class to json and not Geojson.Is there any tool in arcgis desktop to convert shapefile or featureclass into json.


texts: How to get indoor directions by using Google API

texts:
I have a problem using Google Direction API to get the route between two points(One point is outside the building and the other is inside the building).

When I use API to get the directions, only the routes outside building are returned as JSON data. But I want to show the indoor route for my application.

Is it possible to get the indoor directions using this Directions API or other API?


texts: Splitting LineString into Multiple Substrings using PostGIS?

texts:
I am new to PostGIS.

I have a set of (thousands) of Linestrings and associated multiple points along the line, expressed as fractions along each linestring using ST_LineLocatePoint.
I am looking for a wrapper of sorts to take in a Linestring and a list of fractions, and return the geometry collection of all of the created Substrings.

Example input:

 LINESTRING=example_line, FRACTIONS=(0.25,0.835,0.9)

Example Output:

GEOMETRYCOLLECTION(
    ST_Line_Substring(example_line, 0, 0.25),
    ST_Line_Substring(example_line, 0.25, 0.835),
    ST_Line_Substring(example_line, 0.835, 0.9),
    ST_Line_Substring(example_line, 0.9, 1)
)

Does this kind of thing exist in PostGIS?


texts: Dividing a buffer into segments based on path and distance

texts:
I have a buffer layer which I generated from a polyline layer and then dissolved. I then rasterized the dissolved buffer layer, thinned it based on this link and generated a new (more approximate) polyline layer which occurs at more or less the centerline of the dissolved buffer layer.

I also have a background raster containing values which I am interested in averaging later on. In summary my session looks something like this:

I would like to now find a way to segment that buffer layer based on a particular distance. This is the reason I thought that having a centerline polyline layer would be beneficial. In MS paint I have drawn an example of what I am aiming for:

I'd like to make the segments as perpendicular as possible to the edges of the buffer zone. I realize there is also the problem of how the segments are calculated when the road branches off as shown below. To this I would apply the algorithm employed by v.split.length, if possible.

How do I segment the buffer layer based on a regular centerline length of the dissolved polyline layer?


texts: What does the values in column clazz (osm2po) mean?

texts: I have a given table in postgreSQL with OSM data which is used for routing purposes in Java with PGRouting. There exist two columns which I don't understand. They belong to osm2po. They are labeled with "flags" and "clazz". I know, that flags describe the kind of road. It means, whether it is driveable by car, bike,... The problem now is, that I do not understand the meanings of the values in the column clazz. I just know that they also describe the road type. Can anybody explain, what the values like 12, 15, 31,... mean? Does a documentation exist where the meaning of the values is written down?


texts: Stand-alone Server Object Extension (SOE)

texts:
I would like to setup a REST service through ArcGIS Server.

In this case I do not need a map service, but it seems like you can only deploy an SOE as an extension to another service.  It seems weird to attach the SOE to a map service that it is not related to.

I would prefer to create the service as an SOE since it uses ArcObjects and so that it can be managed and deployed through the AGS management console.

I'm working with ArcGIS Server 10.1 on Windows with IIS.