texts: Is it possible to delete geometry (create NULL) but retain the feature row?

texts:
I'm trying to perform the reverse action to "display x/y". I have a point layer displayed on the map canvas, and I want to delete a point's geometry, as in remove it from the map, but retain the feature row and attributes in the layer's table.

I don't want to use any filters or symbology manipulations, just delete its physical geometry

Edit:
Both methods worked, the in-tool expression is slightly more elegant IMO.
the issues raised in the comments abou the validity of such a thing are correct, I should have said that this is a temporary solution during the process of updating point locations and working with a non-GIS person, so I had to, at least temporarily, retain all the features even if they aren't presented on a map.


texts: How to display a continuous diagonal line with an SVG externalgraphic in Geoserver?

texts:
I need to symbolize a polygon layer with geoserver using an SVG file as an external graphic. The polygon should be filled with diagonal line. I have created the SVG in Inkscape using a 400px x 400px SVG Graphic. Geoserver renders the polyons as expected, but the SVG pattern has a gap between them. How can I remove this unsightly gap? This Problem does not occur when served with DEEGREE, so I am assuming it is a geoserver bug.

Link to Layer Hier

Here is the SVG:

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

<svg width="400" height="400" xmlns="http://www.w3.org/2000/svg">;

<polygon points="0 0  50 0 0 50" stroke="#FF0000" stroke-width="1px"
fill="#FF0000" />

<polygon points="150 0 250 0 0 250 0 150" stroke="#FF0000" stroke-
width="1px" fill="#FF0000" />

<polygon points="0 400 0 350 350 0 400 0 400 50 50 400" stroke="#FF0000"
stroke-width="1px" fill="#FF0000" />

<polygon points="400 150 400 250 250 400 150 400" stroke="#FF0000" stroke-
width="1px" fill="#FF0000" />

<polygon points="400 350 400 400 350 400" stroke="#FF0000" stroke-
width="1px" fill="#FF0000" />

</svg>

Here is the SLD Symbolizer.

                <sld:PolygonSymbolizer>
                <sld:Fill>
                    <sld:GraphicFill>
                        <sld:Graphic>
                            <sld:ExternalGraphic>
                                <sld:OnlineResource
                                    xmlns:xlink="http://www.w3.org/1999/xlink"; xlink:type="simple" xlink:href="ROK_gewerbe.svg"/>
                                    <sld:Format>image/svg</sld:Format>
                                </sld:ExternalGraphic>
                                <sld:Size>20</sld:Size>
                                <sld:Rotation>0</sld:Rotation>
                            </sld:Graphic>
                        </sld:GraphicFill>
                    </sld:Fill>
                    <sld:Stroke>
                        <sld:CssParameter name="stroke">#666666</sld:CssParameter>
                        <sld:CssParameter name="stroke-width">1</sld:CssParameter>
                        <sld:CssParameter name="stroke-opacity">1</sld:CssParameter>
                    </sld:Stroke>
                </sld:PolygonSymbolizer>


texts: Cutting polygon using line - cutter, cut() - using ArcPy?

texts:
What I'd like to do, using ArcPy: use a line to split a bunch of polygons into separate parts.
Here is what the inputs could look like:

In this case the line should split each polygon into two separate polygons, creating 6 separate polygons in the end.

This question has been asked many times, and there have been a lot of answers provided, but none of them seem to mention the polygon's cut method, perhaps because it seems to only have been introduced in 10.2.

Here is what the ArcGIS Help shows (scroll to cut (cutter)):

This seems to be exactly what I need, but I cannot entirely figure out how to implement it.

From the few examples I was able to find, it seems like I cannot use the method directly on a shapefile that I am referencing, but I need to access the underlying geometry/vertices, so here is what I did (please correct me if there is an entirely different way of doing it):

read shapefile's vertices using a cursor
push each vertex into an array
use that array to create a polygon/line

The next step would be to to use the cut() method, as shown in the screenshot above, to actually split the polygons, but I cannot get it to work.

What follows is my code (only the polygon-part uses comments as the polyline is just repetition):

import arcpy

polygon = r"D:\Data\poly.shp"
line = r"D:\Data\line.shp"

# ------------------- POLYGON -------------------
# first read the geometry
# i.e.: read all of the vertex coordinates, and push them into an array
dataset_vertices = []

# loop through each feature
for row in arcpy.da.SearchCursor(polygon, ["OID@", "SHAPE@"]):
    print("Here comes feature {0}:".format(row[0]))
    partnum = 0
    # Step through each part of the feature
    for part in row[1]:
        # Print the part number
        print("Part {0}:".format(partnum))

        tempPart = [] # push each vertex into an array
        # Step through each vertex in the feature
        for pnt in part:
            if pnt:
                # Print x,y coordinates of current point
                print("{0}, {1}".format(pnt.X, pnt.Y))
                tempPart.append([pnt.X,pnt.Y])
            else:
                # If pnt is None, this represents an interior ring
                print("Interior Ring:")
        # push part-array into entire array
        dataset_vertices.append(tempPart)
        partnum += 1

# A list that will hold each the Polygon objects
polygon_features = []

for feature in dataset_vertices:
    # Create a Polygon object based on the array of points
    # Append to the list of Polygon objects
    polygon_features.append(
        arcpy.Polygon(
            arcpy.Array([arcpy.Point(*coords) for coords in feature])))

# ------------------- POLYLINE -------------------
line_vertices = []

for row in arcpy.da.SearchCursor(line, ["OID@", "SHAPE@"]):
    print("Here comes feature {0}:".format(row[0]))

    partnum = 0
    for part in row[1]:
        print("Part {0}:".format(partnum))

        tempPart = []
        for pnt in part:
            if pnt:
                print("{0}, {1}".format(pnt.X, pnt.Y))
                tempPart.append([pnt.X,pnt.Y])
            else:
                print("Interior Ring:")
        line_vertices.append(tempPart)
        partnum += 1

line_features = []

for feature in dataset_vertices:
    line_features.append(
        arcpy.Polygon(
            arcpy.Array([arcpy.Point(*coords) for coords in feature])))

After running this code I can use methods, such as polygon_features[0].boundary(), but I am not sure how to use cut().


texts: Deciding what interpolation method to use for resampling raster data?

texts: Are there hard and fast rules about which interpolation methods are suited to each kind of raster data?


texts: How to uncheck toggle button in mapbasic?

texts:
I want to check un check togglebutton but I don't know how to alter using mapbasic. Is there any command to alter it?  I could not find it in MapInfo user guide. I have command:

Alter Button

But how to use I don't know how to provide Id to buttons?


texts:
I am trying to run this from QGIS python console to extract the channel network from the DEM, based on the catchment area:

processing.runalg('saga:channelnetwork', 'C:/Data/filled.tif', None, 'C:/Data/catchment_area.tif', 2, 10000000, None, 10, None, 10, 'C:/Data/channel_network.tif', 'C:/Data/africa/channel_route.tif', 'C:/Data/channel_network.shp')

The output is empty with nan values.

This is how I have derived the catchment area:

import processing
processing.runalg('saga:fillsinks', "C:/Data/clipped.tif", 0.010000, "C:/Data/filled.tif")
processing.runalg('saga:flowaccumulationtopdown', 'C:/Data/filled.tif', 0, 'C:/Data/catchment_area.tif')

I am using QGIS 2.18.11 with SAGA 2.3.2


texts: Error when exporting SOAP Server Object Extension from Eclipse plugin

texts:
I am new to creating SOEs. Limited experience with ArcObjects. We are upgrading from ArcGIS server 10.05 to 10.3 and I am responsible for wrapping all of our code (Java) in SOEs for the upgrade since the code currently uses a DCOM connection. I have found very little documentation on the specifics of writing SOEs but have installed the Eclipse plugin (Ver 10.3.1.4959) in Rational Application Developer (RAD 9.6) as instructed to be able to create and export the SOEs. The documentation that I have been following is ESRI’s ArcObjects Help for Java – Developing Extensions.

I am getting errors when trying to export a SOAP SOE using the plugin. It appears to export everything else but fails when it is exporting the actual soe file. Error: “javac” terminated with errors.

The one resolution for this that I found is setting your PATH variable. I tried that, it doesn’t resolve my issue.  Is there anything else I can do to fix this?

NOTE: I have also asked this question on GeoNet community and have not received any response.


texts:
I have some code in loop:

var vector =   new ol.layer.Vector({
                                title: name,
                                style: styleFunction
                            });
                             function getFeatures(geoJson) {
                                // If response is valid
                                var format = new ol.format.GeoJSON({
                                    featureProjection: "EPSG:3857"
                                });
                                vector.setSource(
                                        new ol.source.Vector({
                                            format: new ol.format.GeoJSON(),
                                            features: format.readFeatures(geoJson),

                                        }));
                            }

                            $.ajax({
                                type: "GET",
                                url: urlAjax,
                                dataType: 'jsonp',
                                jsonpCallback: 'getJson',
                                success:getFeatures
                            });

                            overlays.push(vector)

After push vector to overlay array I add that array of vectors to map:

 var map = new ol.Map({
        target: 'map',
        layers: [
            new ol.layer.Group({
                title: 'Base maps',
                layers: basemaps
            }),
            new ol.layer.Group({
                title: 'Overlays',
                layers: overlays
            })

        ],
        view: new ol.View({
            center: ol.proj.fromLonLat([31.4, 48.9]),
            zoom: 6
        })
    });

And create switcher controls... And all vectors are reset sources last layer and only last layer visible in map, others don`t work...why?


texts: In ArcMap graticule ticks on opposite ends of the dataframe do not align

texts: I'm adding tick marks and labels to my dataframe but the ticks on the top do not align with the corresponding ticks on the bottom.  The same is true with left and right ticks.  See image below for visual aid. I want to create a map with three adjacent dataframes but it is difficult to label efficiently if the ticks are not aligned.  My best guess is this is because of the coordinate system and the fact that I'm attempting to display a curved surface in a flat box.  Is this my problem or is there another solution?  I'm using ArcMap 10.3.1 and the coordinate system is coordinate system is NAD 1983 UTM Zone 18N.


texts:
I'm looking for a metric that gets at average distance to a given feature. For example, for a point in the water, how far, on average, is the land or to what extent is the point surrounded by land.

Scenario: two points are the same distance from shore, but one is adjacent to the tip of a narrow peninsula of land while the other is at the bottom of a narrow bay. Clearly one point has more land around it than the other. I'm looking for a way to quantify this.

I'm not concerned with the specific software tool. Formulas or names of metric would be useful to get me looking in the right place in the literature.


texts: Projection to use with satellite image from Google Maps / Earth

texts:
I'm very new to anything concerning GIS so please forgive me if this question has been asked before, but I assure you, I've been trying to find an answer for quite some time now. Of course it is possible, that the question itself is dumb or that I've overlooked the answer because of my lack of expertise in this field.

I built a programme for a 3D Engine that reads latitude and longitude values from an openstreetmap (osm) file and projects these onto a x/y plane. With the projected corner points of buildings I then create rudimentary 3D buildings by just extruding their ground areas. I used two projections so far, the Mercator and the UTM projection.

Now I'd like to "download" satellite images (or rather make screenshots from Google Maps / Earth) and use the as an underlay, so that the constructed buildings sit on top of the image in the right places and the right scale. I'd like to be able to align the buildings with the buildings on the satellite image by doing no more than the following three actions (since I have to do it manually):

Moving the image in the x/y plane
Rotating it around the z axis
Scaling it (maybe independently) in x and y

I cannot, however, distort the image in any other way. I know that with each projection you have a tradeoff, as there is no projection that is equidistant, isogonal and where the area proportions are correct all at the same time.

With the two projections I've used so far it wasn't possible to align the image with the buildings, since as soon as I made it fit in one place it didn't fit in other places any more.

So my question is, which projection I should use and which properties the projections needs for my purpose. As I see it, a satellite image (if what you get from Google Maps / Earth even is from a satellite) is a projection itself, so is there are projection that is similar to what a satellite sees when orbiting the earth?! Are there any others sources for satellite images that would work better?

Thanks in advance!

P.S. The area covered by the osm file is rather small (around 3km x 5km).


texts: RESTful Geoprocessing Service Request returns empty Result object without errors

texts:
This is a bit difficult for anyone here to test as our ArcGIS Server 10.2 instance is behind a firewall.  I'll try to present all the necessary information, but feel free to ask for more and I'll adjust the question accordingly.

I have a Geoprocessing Service hosted on ArcGIS Server.  The service is a very simple Python script (below) that takes in a string of Hydrologic Unit Codes (in this case sub-watersheds [HUC8 level]).  The script takes this string, splits it and compares against a table of species occurrence by watershed.  This data is compiled into a Python dictionary and returned as JSON with the json module.

I'm accessing this service via an ajax request from a web mapping app, but for simplicity's sake I'll leave the JavaScript side out unless it's really necessary.

Python

import arcpy, json

# files
hucs = r'S:\ES\ES 2007 Efiles\PROP-500 Maps and Photography\PROP-540 Geographical Information System (GIS) Data files\MegaPetition.gdb\HUCs'
species = r'S:\ES\ES 2007 Efiles\PROP-500 Maps and Photography\PROP-540 Geographical Information System (GIS) Data files\MegaPetition.gdb\CandidatesHUC'
data = {}
Watersheds = arcpy.GetParameterAsText(0)
Watersheds = Watersheds.split(",")

cursor = arcpy.SearchCursor(hucs)

for row in cursor:
    for watershed in Watersheds:
        if row.HUC_CODE == watershed:
            currentWatershed = {}
            speciesList = []
            arcpy.AddMessage("Working on " + row.HUC_NAME + "...")

            whereClause = "HUC8_CD = '" + watershed + "'"
            speciesCursor = arcpy.SearchCursor(species, whereClause)
            arcpy.AddMessage("    Contains " + str(row.NumPetitionedSp) + " At-Risk Species.")

            for animal in speciesCursor:
                speciesList.append(animal.GCOMNAME)
            currentWatershed = {'attributes': [watershed, row.NumPetitionedSp],
                                'species': speciesList}
            data[row.HUC_NAME] = currentWatershed

print json.dumps(data)
del row, cursor

Perhaps I don't understand how the script gets the string of watersheds from the request so that I can use it in my script.  With a Script Tool I use arcpy.GetParameterAsText(0), but that doesn't seem to apply here.

I've tried to execute this script through the ArcGIS Server GP Service Page, but I still get an empty result:

Generates the following URL:
http://server-address/arcgis/rest/services/NatureServe/AtRiskGPService/GPServer/NatureServe%20At%20Risk%20GP%20Service/execute?Watersheds=06040003%2C06040002%2C06030003%2C06030002%2C06030005&env%3AoutSR=&env%3AprocessSR=&returnZ=false&returnM=false&f=pjson

With the following Response:

Update!

I changed the messaging option for the service from Error to Info, so I can see that the script is taking in the watershed parameter and processing it.  It looks like it's just not sending the response back as part of the Result object.

The Result object is created by running a GP tool (ex. result = arcpy.GetCount_management(in_table)).  Is there any way I can add my own data to the Result object if I'm not actually running any GP tools?  If not I guess I need to send the JSON through a message instead.


texts: A loop to select all active layers in a project

texts:
I am trying to export all active layers in a QGIS project to a specified folder, below is what I get atm

Outputdir='C:/Users/DN/Downloads/Test/'

layer = iface.activeLayer()

QgsVectorFileWriter.writeAsVectorFormat( layer, Outputdir + layer.name() + ".shp", "utf-8", layer.crs(), "ESRI Shapefile", 1)

It only works when I click on the layer that I want to export, I am wondering if it's possible to write a loop to loop through all active layers in the project and export them.


texts: Verifying that ArcInfo Workstation is Installed?

texts: What is the best way to verify that a machine has ArcInfo Workstation installed? A code sample in VBA, VB or C# would be great!


texts: Identifying duplicates between two different feature classes

texts:
I am trying to identify duplicate records between two feature classes in a file gdb (Arc 10.3). I first create a set of unique IDs present in FC_1, then use an if statement to identify whether the IDs in FC_2 are present in FC_1.

checkfield = "gisID"
unique_values = set(row[0] for row in arcpy.da.SearchCursor(FC_1, checkfield))

with arcpy.da.UpdateCursor (FC_2, [checkfield, "Type"]) as uCurs:
        for uRow in uCurs:
            print (uRow[0])
            if uRow [0] in unique_values:
                uRow [1] = "duplicate"
                uCurs.updateRow(uRow)
            else:
                uRow [1] = "unique"
                uCurs.updateRow(uRow)

    # check if uRow[0] is in the set
    if 200683 in unique_values:
        print ("200683 is in the set")

I've purposefully set up the feature classes such that the 2 records in FC_2 have IDs already present in FC_1. The last 3 lines in the script are there to check this.

However, FC_2 "Type" field only updates with "unique", when these are clearly duplicates. Any advice on how to correct this?


texts: How do I export neighbourhoods in wikimapia to KML file

texts:
Possible Duplicate:
  How do I export from Wikimapia to KML

I am trying to export the neighbourhoods in Wikimapia to KML format. Is this possible? These are all the predefined neighbourhoods that's gets highlighted on the mouse over.


texts:
I am attempting to convert vector files to raster format in QGIS using the Rasterize tool and am receiving the error:

  ERROR 1: Wrong value for -tr parameter.

I have attached an image of the input and errors I am getting. And I am using QGIS v 3.2.3.


texts: Displaying SHOM WMS in Folium / Leaflet map

texts:
I'm trying to display this WMS in Leaflet / Folium: http://services.data.shom.fr/INSPIRE/wms/r

I use the following code:

import folium

center=[49.9765229931438, 1.1902348134956244]

def inline_map(m, width=800, height=500):
    """Takes a folium instance and embed HTML."""
    m._build_map()
    srcdoc = m.HTML.replace('"', '&quot;')
    embed = HTML('<iframe srcdoc="{}" '
                 'style="width: {}px; height: {}px; '
                 'border: none"></iframe>'.format(srcdoc, width, height))
    return embed

width, height = 1400, 800

#tileset = r'http://server.arcgisonline.com/ArcGIS/rest/services/Ocean_Basemap/MapServer/tile/{z}/{y}/{x}'

tileset = r'http://services.data.shom.fr/INSPIRE/wms/r/{z}/{y}/{x}'

map = folium.Map(location=center, zoom_start=10,
                    #tiles='Stamen Terrain',
                    tiles = tileset,
                    #tiles='OpenStreetMap',
                    width=width,
                    height=height,
                    attr="MyDataAttribution")

map.save("mapTEST.html")

When I open mapTEST.html with my brother it does create a leaflet map object, but no tile is displayed, when I use another WMS, it works perfectly. Did I do something wrong with the URL?


texts: How to get the number of polygons of a geometry in PostgreSQL

texts:
How can I get the number of polygons that form a geometry in PostgreSQL?

For example U.S map has a lot of islands in addition to its main lands. I want to count these islands plus the main lands.


texts: How can I show photos in popup from a ByteA (postgis) in QGIS (or Google Earth/map)

texts: I currently have a postgis database with a table containing point geometry of fire Hydrant. This table contains a raster_blob field which contains the Fire hydrant picture taken by field workers.  Would there be a way to read the picture from the PostGIS ByteArray so it can be shown in a label or pop-up using QGIS (or Google Earth/map)?


texts: What is the equation needed for forward and inverse transformation of Geostationary Satellite projection.?. How to convert it to Geodetic and vice versa. Geostationary satellite projection is mentioned in Proj4 as here and in wiki


texts: equivalent of sp package's point in polygon / overing using sf

texts:
I am migrating code from sp package to the newer sf package. My previous code I had a polygon SpatialDataFrame (censimentoMap) and a SpatialPointDataFrame (indirizzi.sp) and I got the polygon cell id ("Cell110") for each point laying within with the instruction below:

points.data <- over(indirizzi.sp, censimentoMap[,"Cell110"])

Actually I created two sf objects:

shape_sf <- st_read(dsn = shape_dsn)
shape_sf <- st_transform(x=shape_sf, crs=crs_string)
and
indirizzi_sf = st_as_sf(df, coords = c("lng", "lat"), crs = crs_string)

And I am looking for the sf equivalent of the above instruction...
Migth it be:

ids<-sapply(st_intersects(x=indirizzi_sf,y=shshape_sfpeCrif), function(z) if (length(z)==0) NA_integer_ else z[1])
cell_ids <- shape_sf[ids,"Cell110"]


texts:
In pursuit of replicability (as well as requirements of funding agency), I am trying to produce all of the maps for a current project in R.

For single attribute maps, so far so good! I am using ggplot2 and enjoying all of the customization options.

However, I'm running into problems as to how I can make multiple attribute maps. In GIS terms, I'm interested in having one graduated color layer and one proportional symbol layer. Here is an example of the current code:

mapPlot <- ggplot(data) + theme_void() + theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank() ,axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), legend.position="bottom", legend.key.width = unit(1, "cm"), legend.text = element_text(size = 8))

pdf("path.pdf",width=7,height=5)

map <- mapPlot + geom_polygon(data=data, aes(x = long, y = lat, group = group, fill = data$var1), lwd = 1/10, color="lavenderblush4") +
  scale_fill_continuous(name = "Average", low="white", high="dodgerblue4", guide="colorbar", limits=c(0.0,1.0), labels=c("6", "661", "1,106", "2,037", "11,846")) +
  theme(legend.text = element_text(size=10))

map <- map + coord_map(project="conic", lat0 = 30)
plot(map)

This is producing a nice graduated colors map, but I'm unsure how to combine it with an additional layer (and thus additional legend).

I would like to something (relatively) straightforward such as asked here: How do I plot points as graduated/proportional circles in R?

EDIT: Information about my data

-Currently using commuter zone areal units

-Do not have centroids for these units (yet - would they help?)

EDIT 2: Found this link (https://datacarpentry.org/r-raster-vector-geospatial/08-vector-plot-shapefiles-custom-legend/), going to update my results

Edit 3: No luck via link.


texts:
I have a set of points which represents a series of telegraph poles. I want to link them to simulate the wire running along them. The code I've been using is shown below.
SELECT DISTINCT ON (p1.id)
       p1.id AS p1id,
       p2.id AS p2id,
       ST_ShortestLine(p1.geom,p2.geom) As snapped_point
FROM public.poles AS p1
INNER JOIN public.polescopy AS p2 ON ST_DWithin(p1.geom,p2.geom,100.0)
WHERE p1.id <> p2.id
ORDER BY p1.id

A problem occurs when the next pole in line is further away than the previous pole. This causes the nearest neighbour search to double back to the previous pole and ends the line prematurely.
Is there a way of tagging the previously used pole so that the nearest neighbour search ignores it?


texts: Emulating ArcGIS Highest Position tool in QGIS?

texts:
I am familiar with the tool Highest Position (Spatial Analyst) in Arcgis, but I am trying to replicate this tools function is QGIS. Is anyone aware if this tool exist in QGIS (possibly under a different name)? I have tried looking but no luck yet. Below is a description of what the tool does in Arcgis.

So the way the Highest Position Tool works in Arcgis is as follows:

The order that you added your raster is important for the output:
When I add InRas1 to the tool, it is my first input for my stack of rasters. The tool assigns this raster the value of 1
When I add InRas2 to the tool, it is my 2nd input for my stack of rasters. The tool assigns this raster the value of 2
When I add InRas3 to the tool, it is my third input for my stack of rasters. The tool assigns this raster the value of 3.
When the Highest Position tool is run, it is looking at over laying pixels to determine the highest value. In this example I will use the bottom left pixel of inputs InRas1, InRas2 and InRas3. From the illustration, InRas1 has the highest value of 4 and so in the output raster (outRas), the tool assigns that pixel a value of 1 (as this was the first raster in putted into my raster stack. See point 2 above)


texts: Python script/Model to populate cross street names at intersections

texts:
Has anyone created a Python version/model of this VBA script from ArcScripts which populates the to and from street names at street intersections?  A logical workflow is outlined quite nicely in this ESRI thread, but it requires ArcInfo, and I would like to know if there is another way to do it.

The end result would be, for example, if there are lines intersecting at First Street and Third Street, and First Street and Main Street, the record for First Street would have to = Third Street and from = Main Street (the order does not matter).

Using Select By Location to select connected roads seems like it would work, and could even possibly replace the need to create points in the aforementioned workflow (looping over each street, writing to separate fields the streets it touches based on the output of the select by location, then extracting the two street names which are different). Does anyone have any other ideas?


texts: Different Behaviour of GDAL.RasterIo with Python and C++

texts:
I was using GDAL C++ API and Python API to open a tif dataset and read it into memory, but then I got some annoying issue, when I use Python,all works fine.

import gdal
ds = gdal.Open('/*.tif')
band = ds.GetRasterBand(1)
nXsize = band.XSize
data = band.ReadRaster(xoff=0, yoff=0, xsize=nXsize, ysize =1, buf_xsize=nXsize, buf_ysize=1)

And the out put is right, nothing happend.

data: b'\xe6\x02\xce\x02\xc7\x02\x8e\x02B\x02\xfc\x01\x12\x02\xf5\x01\xd2\x01\xc3\x01\x04\x02\xfc\x01\xc9\x01\x90\x02\xf9\x02\xe4\x02p\x02\x1b\x02\xe8\x01\x02\x02\x00\x0...

Then comes C++ API

#include <iostream>
#include <gdal_priv.h>
#include <cpl_conv.h>

using namespace std;

int main()
{
    const char* pszFile;
    GDALAllRegister();
    pszFile = "/*.tif";
    GDALDataset* poDataset = (GDALDataset*)GDALOpen(pszFile, GA_ReadOnly);
    GDALRasterBand* poBand = poDataset->GetRasterBand(1);
    float* pafScanline;
    int nXSize = poBand->GetXSize();
        pafScanline = (float*)CPLMalloc(sizeof(float) * nXSize);
    poBand->RasterIO(GF_Read, 0, 0, nXSize, 1,
        pafScanline, nXSize, 1, GDT_Float32,
        0, 0);
    system("pause");
    return 0;
}

I comliles the program, nothing happend, then I run it, errors comes out:

ERROR 6: Unhandled type for StripOffset/StripByteCount
ERROR 1: /*.tif, band 1: IReadBlock failed at X offset 0, Y offset 0: Unhandled type for StripOffset/StripByteCount

The xoff,yoff,Xsize,Ysize parameter are same in two API, why C++ gives me error output?The error said itIReadBlock failed at X offset 0, Y offset 0, isn't the xoffset and yoffset start with 0?Hope somebody could answer my question.


texts: Leaflet's shapefile plugin only compatible with Firefox

texts:
I created an HTML page (currently locally) and used Leaflet to create a basic map to which I plan to add several shapefiles. To do this, I use the leaflet.shapefile extension (https://github.com/calvinmetcalf/leaflet.shapefile) and everything is working pretty good so far.

That is, in Firefox. If I test the page in Internet Explorer or Google Chrome, everything except for the shapefile I use for testing shows up.

I uploaded the page to mediafire (http://www.mediafire.com/download/wdib96x6udbnq2h/web.rar), index.html is the main page.

I consider myself pretty seasoned in GIS and Cartography terms, but I'm fairly green when it comes to the web stuff like HTML, so if this not the right question for this forum feel free to point it out : ) .


texts: Problems with labels halo in OpenLayers2

texts:
Draw label as follows:

var styleMarkersLabels = new OpenLayers.Style(
                {
                    strokeWidth: 2,
                    labelYOffset: 14,
                    label: "${label}",
                    fontColor: 'red',
                    fontSize: 16,
                    fontWeight: "bold",
                    labelOutlineColor: "black",
                    labelOutlineWidth: 1
                });
            var vectorPoint = new OpenLayers.Layer.Vector("Узлы (надписи)",
            {
                styleMap: new OpenLayers.StyleMap(
                { "default": styleMarkersLabels,

                "select": { pointRadius: 20}
                })
            });
            map.addLayer(vectorPoint);

            function addPoint(lon, lat, title, ident, layr){
                var ttt = new OpenLayers.LonLat(parseFloat(lon), parseFloat(lat));
                ttt.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));
                for (var k = 0; k < layr.features.length; k++) {
                    if(layr.features[k].attributes.PointId==ident) {
                        layr.features[k].move(ttt);
                        layr.features[k].attributes.label=title;
                        return false;
                    }
                }
                var point0 = new OpenLayers.Geometry.Point(parseFloat(lon), parseFloat(lat));
                point0.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));
                layr.addFeatures(new OpenLayers.Feature.Vector(point0, { label: title, name: title, PointId: ident }));
            }

            var lat2, lon2, title, ident;
            for (l = 0; l < nodesLabels_Count; l++) {
                lat2 = nodesLabels_arr[l]["points"][0]["lat"];
                lon2 = nodesLabels_arr[l]["points"][0]["lon"];
                title = nodesLabels_arr[l]["title"];
                ident = nodesLabels_arr[l]["ident"];
                addPoint(lon2, lat2, title, ident, map.layers[4]); // Here adds the label
            }

But the halo is shifted slightly. Screenshot: http://softez.pp.ua/Image.PNG

How to fix it?


texts: Drawing line or polygon interactively in ArcMap and later saving it as shapefile?

texts:
I am using ArcGIS 10.4 for Desktop.

I need to draw line feature in current mapdocument interactively and later save it for further analysis such as stack profile so that I can do the stack profile generation interactively.


texts:
I would like to delete all rows which are joined to a layer from a csv table. Only the joined records which match (keep only matching records) should be deleted. After the records are deleted, the join should be removed again. I used the following script, but get the error: "cannot update join table".

arcpy.env.workspace = r"D:\skript\New File Geodatabase.gdb\test"

arcpy.MakeFeatureLayer_management("fc","featurelayer")

arcpy.AddJoin_management("featurelayer","OBJECTID",r"D:\skript\Errors.csv","id","KEEP_COMMON")

with arcpy.da.UpdateCursor("featurelayer", "OBJECTID") as cursor:
for row in cursor:
    if row[0] > 0:
        cursor.deleteRow()

arcpy.RemoveJoin_management ("featurelayer", "Errors")

arcpy.CopyFeatures_management("featurelayer","result")

Is this even possible?


texts: Convert From Northing/Eastings to Latitude/Longitude with QGIS

texts:
How can I convert a file with the projection below (NAD83) to one with Latitude and Longitude?

The import file is a SHP file and I would like to export back out a SHP file.

+proj=tmerc +lat_0=0 +lon_0=-114 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

A bit of a beginner with GIS and QGIS but familiar with ArcGIS though I no longer have access to it.


texts: struggling with GeoDjango install: ERROR: could not access file "$libdir/postgis-1.5": No such file

texts:
I'm installing GeoDjango on Mac OSX. I've followed the official Mac install instructions, which all went smoothly, and am now creating a spatial database template for PostGIS.

However, when I try to load the PostGIS SQL routines, I get ERROR:  could not access file "$libdir/postgis-1.5": No such file or directory:

postgres$ psql -d template_postgis -f $POSTGIS_SQL_PATH/postgis.sql
psql:/usr/local/pgsql/share/contrib/postgis-1.5/postgis.sql:59: ERROR:  could not access file "$libdir/postgis-1.5": No such file or directory
<snip>
psql:/usr/local/pgsql/share/contrib/postgis-1.5/postgis.sql:7785: ERROR:  type "geometry" does not exist

What's going wrong, and how can I fix it?

I think I'm running the right version of pg_config, etc:

postgres$ which pg_config
/usr/local/pgsql/bin/pg_config
postgres$ which psql
/usr/local/pgsql/bin/psql
postgres$ pg_config --pkglibdir
/usr/local/pgsql/lib

If I look in /usr/local/pgsql/lib, there is a file called postgis-1.5.so there, so I don't understand why it's not found.

-----UPDATE-------

Does this tell anyone anything?

delirium:~ ap257$ otool -L /usr/local/pgsql/lib/postgis-1.5.so
/usr/local/pgsql/lib/postgis-1.5.so:
        /Library/Frameworks/GEOS.framework/Versions/3/GEOS (compatibility version 4.0.0, current version 4.2.0)
        /Library/Frameworks/PROJ.framework/Versions/4/PROJ (compatibility version 7.0.0, current version 7.6.0)
        /usr/lib/libxml2.2.dylib (compatibility version 10.0.0, current version 10.3.0)
        /usr/lib/libz.1.dylib (compatibility version 1.0.0, current version 1.2.3)
        /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 125.0.0)
        /usr/lib/libicucore.A.dylib (compatibility version 1.0.0, current version 40.0.0)


texts:
I would like to have a map with the Sahel. Namely 'Mali', 'Senegal', 'Guinea', 'Guinea-Bissau', 'Sierra Leone', 'Ivory Coast', 'Mali', 'Niger', 'Chad', 'Sudan', 'Eritrea', 'Ethiopia', 'Benin', 'Burkina Faso', 'Togo', 'Ghana', 'Nigeria', 'Cameroon', 'Gambia', 'Liberia', 'Mauritania' should be included.

The map should be really basic and the country names should be written within the boundaries.

Any idea how to code it with R?


texts: Unable to dynamically locate map boundary coordinates in OSM/Leaflet

texts:
I have been trying for the past hour to take the coordinates of the map edge in a OSM/leaflet application I am making. While I can get the coordinates of the center, and I can get the coordinates of the two poinths the user inputs, the bounds of the map remain unchanging and static on the ones the map loaded. I think it has to do with the map.setView() I call in the start of the code, but I am not sure how to fix that. The code is as follows (the problematic part is actually on the lower part of the code):

var map = L.map('map').setView([37.92,23.69], 13);
var bounds=map.getBounds();
var nodeArray=new Array();
var markerArray = new Array();
var startPoint=null;
var endPoint=null;
var lowerBounds=bounds.getSouthWest();
var upperBounds=bounds.getNorthEast();
var markerCounter=0;

L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © CloudMade',
maxZoom: 18}).addTo(map);

var popup = L.popup();
var road = L.polyline(nodeArray).addTo(map);

map.on('click', onMapClick);
function onMapClick(e) {
    markerCounter=markerArray.length;
    if (markerCounter==2){
        map.removeLayer(markerArray[0]);
        map.removeLayer(markerArray[1]);
        markerArray=new Array();
        endPoint=null;
        startPoint=null;
    }else{
        var marker = new L.Marker(e.latlng, {draggable:false});
        markerArray.push(marker);
        for (i=0;i<=markerArray.length;i++){
            map.addLayer(marker);
        }
        if (markerCounter==1){
            marker.bindPopup("Your destination is: "+(e.latlng.toString()).substr(6,20)).openPopup();
            endPoint=(e.latlng.toString()).substr(6,20);
            markerCounter=2;
        }
        if (markerCounter==0){
            marker.bindPopup("Your starting place is: "+(e.latlng.toString()).substr(6,20)).openPopup();
            startPoint=(e.latlng.toString()).substr(6,20);
            markerCounter=1;
        }

    }
};

function runScript(){
    if (markerCounter==0){
        alert('You have give us no data. We cannot calculate a route');
    }else if (markerCounter==1){
        alert('You have give us only one point. We cannot calculate a route');
    }else{
        console.log('LB:'+lowerBounds+" - "+"UB:"+upperBounds+"/SP:"+startPoint+" - EP:"+endPoint);
        }
}

Does anyone have any idea what is the problem?I am unfortunately entirely new to JS and OSM/Leaflet, so I cannot really see easily problems like those.


texts: Open Data Cube: Aggregation of pixels' value based on shapefile and mask raster

texts:
I have several Sentinel-2 products connected to ODC. My purpose is to use a shapefile with parcels in order to extract the average value of each parcel (creation of object based dataset). In addition, there is a raster (SLC) which shows the probability of cloud on each pixel. Therefore I do want to mask and leave out all the pixels mentioned as clouds.
How can I make the aggregation of the pixels per parcel and at the same time the cloud masking based on SLC (Until now this process takes place with gdal and pure Python and I need to implement this via ODC)?


texts: Remove notification when shp selected automatically in startup.py

texts:
I have a script from List index out of range when put in startup.py on QGIS :

from qgis.core import *
layer =  QgsVectorLayer('C:/Users/path_to_shapefile/vector_layer.shp', 'pl' , "ogr")
QgsMapLayerRegistry.instance().addMapLayer(layer)
vl = QgsMapLayerRegistry.instance().mapLayersByName('pl')[0]
iface.setActiveLayer(vl)

I want to put the script in startup.py. So when load QGIS, 'pl.shp' layer will be selected.

With this script, there is notification to save project or not. I have read to remove notification in Remove saving notification when use python QGIS at startup. Then I add the script in Remove saving notification when use python QGIS at startup. But, there is still notification. So, how to combine two script without saving notification?


texts:
I run the code below, that I would hope to return the parameters from projection file (.prj). However it doesn't work as I would expected. I realized that the it isn't working because the some features of .prj file. I tried with another .prj file and the code has worked well.
What is the problem to my code?

import ogr

# use Shapefile driver
driver = ogr.GetDriverByName("ESRI Shapefile")
# reference Shapefile
shp = "/home/dogosousa/Downloads/lc8_ba.shp"
# open the file
ds = driver.Open(shp, 0)
# reference the only layer in a Shapefile
lyr = ds.GetLayer(0)

print(lyr.GetSpatialRef())

# projected coordinate system
proj_string = lyr.GetSpatialRef().GetAttrValue("PROJCS", 0)
# geographic coordinate system
geog_string = lyr.GetSpatialRef().GetAttrValue("GEOGCS", 0)
# EPSG Code if available
epsg = lyr.GetSpatialRef().GetAttrValue("AUTHORITY", 1)
# datum
datum = lyr.GetSpatialRef().GetAttrValue("DATUM", 0)

print(proj_string, geog_string, epsg,datum)

The .prj file that doesn't work with the code:

PROJCS["PROJCS["WGS 84 / UTM zone 24N",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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32624"]]",UNIT["Meter",1]]

The .prj file that works with the code:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]


texts: Calculating coordinates of square x miles from center point?

texts:
I am trying to create a x mile square (or circle) around a central point, where all sides of the square would be x miles from the centre. I need the 4 corner coordinates.

It is scrambling my brain trying to get my head round it?
I can work out the distance between two points using the haversine formula but maths is seriously not my strong point and I don't understand sin,cos etc.. and trying sort this out has lost me!

I have come across Calculating Latitude/Longitude X miles from point? but I just don't get it!

Would anyone be kind enough to explain how I do this in apples and pears terms?

To explain exactly what I trying to do;

I have a website, where users can search for buildings in a specific area. They will enter a town or place (which I will know the lat long of) and they search within a specific radius of say 10 miles of the place.

I need to find the min/max lat and longs of the 10mile radius so I can query my database using a where clause similar to:

Where buildingLat <= maxLat
  and buildingLat <= minLat
  and buildingLong >= minLong
   or buildingLong >= maxLong

I need some kind of formula!

My coordinates are in decimal degrees


texts:
I have downloaded the raster data for Europe from https://land.copernicus.eu/imagery-in-situ/eu-dem (EU-DEM v1.1).
To these data I apply the following orders:

gdalwarp -s_srs EPSG: 3035 -t_srs EPSG: 3785 -r bilinear dc.tif dc-3785.tif
gdaldem hillshade -co compress=lzw dc-3785.tif dc-hillshade-3785.tif
gdaldem slope dc-3785.tif dc-slope-3785.tif
gdaldem color-relief -co compress = lzw dc-slope-3785.tif slope-ramp.txt dc-slopeshade-3785.tif

to get the raster layer, but when inserting the hillshade-3785.tif and slopeshade-3785.tif layer into my map, the map takes a long time to load the tiles and sometimes doesn't even load them.
How can I make the downloaded raster data not so heavy so that the map loads fast and doesn't slow down?


texts: Publishing Geoprocessing Service with ArcGIS Server?

texts:
I'm modeling this ESRI tutorial for creating geoprocessing service with ArcGIS Server and I'm getting this error when trying to publish the service:

  Service was published but cannot be started due to: Server object
  instance creation failed on all SOC machines.  Server Object instance
  creation failed on machine MyServer.  Cannot load tool: BufferAndClip:
  Geoprocessing service initialization failed.

I'm using ArcGIS Desktop 10 against ArcGIS Server 9.3.  My mxd that has the model in it has been saved out as version 9.3.  Also my geodatabase is a 9.3 version.  The model works fine within the mxd.  The mxd and geodatabases and lyr files are located all in the same folder on the server (e.g. \\MyServer\ArcGIS Resources\Test\Test.mxd).

I have published map services before to this server with no issues, why I'm I getting this error?


texts:
Does anyone know how to convert a .lyr to .shp format in ArcMap 10.2?

I had a look at different forums on the internet and most people suggest that it works by right clicking on the layer in the TOC and export the data as shp, but apparently this is not true because the file extension is *. lyr .

The layer I want to export is a web map.

Does anyone have any suggestions?


texts: How to display a map from MapServer in OpenLayers

texts:
I'm working on a project with MapServer and I need to display my generated map in OpenLayers, but I'm not understanding how to do this. I've searched but none of the results helped.

I'm using OpenLayers 4.4

Either OpenLayers documentations helped.

I've tried these two exemples that I found on OpenLayers page but none of them worked.

http://openlayers.org/en/latest/examples/arcgis-image.html?q=mapserver

http://openlayers.org/en/latest/examples/arcgis-tiled.html?q=mapserver

I've changed the url to use my source but didn't work.


texts: Libraries to navigate and choose path based on congestion level?

texts:
I'm trying to develop a system to get best path to travel user based on road traffic (congestion level). Functionality of the system should be provide least congested path for traveler given destination from his/her current location.

Note - I'm using OSGeo live 6.5 , Postgresql , MapServer wms services etc. for now.

There are routing libraries like pgroute to get shortest path. Is there such libraries/tool that we can use for above stated development?


texts: Labeling text with special characters in ArcMap?

texts:
I am trying to display some bird data in ArcMap in the format/notation that it is recorded in the field.  The notation is made up of the one/two letter species code followed by a description of the activity underneath (see image below).

I do not want to create symbology because I have 500+ species and six possible activities, so I have decided to use labelling in ArcMap rather than symbology.  I am able to present data with a circle and underlined but I am struggling with displaying the data that has a double underline or asterix underneath.  Therefore, I was wondering if anybody would be able to offer me any help - how do I include these characters in a label?  I am running ArcMap 9.3.


texts: Why do R scripts not show up amongst the Providers in QGIS3.18?

texts: I have installed the Resource Sharing plugin (3.18) and have installed the QGIS R script collection as instructed elsewhere. Why do the R scripts not now show up among the Providers or in the Processing toolbox?


texts:
Is there a simple way to generate labels for a multiband Raster in QGIS like it is implemented as Default in ArcGIS?
I want to write a Variable Name of each Band in the legend
ArcGIS example:


texts: Are there another way/tool as GDAL/ERDAS-API to convert a raster to ECW?

texts:
I am looking for a tool to convert raster (TIF, JPG) to ECW. I found this two year old question/answer How to convert orthophotomap to .ecw file? The answer is to use GDAL with an ERDAS SDK. I checked the current free versions of ERDAS SDK: The free version of ERDAS SDK is still read-only.

Is there meanwhile another way/tool as GDAL/ERDAS-API to convert a raster to ECW?

(I prefer a free program. If it is not free, information about the approximate cost is helpful.)


texts: Transform GML to Geometry (PostGIS) in GeoKettle

texts:
How can I transform GML (saved in XML file) to Geometry type?

I've noticed function 'Geometry from GML' in Calculator but it needs string value with GML. Is there any way to extract it from XML file? 'Get data from XML' component can give me values stored in node but not whole node.


texts: How to apply a cloud mask on a MODIS collection in Earth Engine

texts:
I am trying to work with the MODIS dataset from google earth engine (MCD43A4.006 MODIS Nadir BRDF-Adjusted Reflectance Daily 500m) and I want to apply a cloud mask but I have no clue how to do it. I also want to keep the images with 30% or less of clouds (masked pixels). I have done this manipulation before with a landsat collection but I am not able to do it with MODIS.
Did anyone know how to do it ?


texts: Get filepath then append a feature class to that feature class Python Addins

texts:
I'm trying to use tool bar buttons (python addins) to add a filepath and then append another feature class (known to have identical schema, no test needed) to that filepath. My syntax isn't working for some reason.

import arcpy, os
from arcpy import env
import pythonaddins
import re
sdeconn1 = r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\Connection to ffxsde.sde"

class WORKSPACEPATH(object):
    """Implementation for addin_addin.WORKSPACEPATH (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        layer_files = pythonaddins.OpenDialog('Select Layers', True, r'C:\GISData', 'Add')
        mxd = arcpy.mapping.MapDocument('current')
        df = pythonaddins.GetSelectedTOCLayerOrDataFrame()
        list22 = []
        if not isinstance(df, arcpy.mapping.Layer):
            for layer_file in layer_files:
                layer = arcpy.mapping.Layer(layer_file)
                list22.append(layer_file)
                pythonaddins.MessageBox(str(list22), 'INFO', 0)
        else:
            pythonaddins.MessageBox('Select a data frame', 'INFO', 0)

class APPEND(object):
    """Implementation for addin_addin.APPEND (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        arcpy.Append_management(r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb\STW_ST_SIGNS", list22, "NO_TEST","","")


texts: How do I run Panoply (from NASA GISS) on Windows 10?

texts:
I am unable to run Panoply (from NASA GISS) to view NetCDF files. When I try to run it, it tells me that I need Java Runtime Environment 9, but the latest version of JAVA from the Oracle website is 8.271. How do I resolve this? I have found no documentation for the same. I want to use it to visualize MODIS Leaf Area Index data.
I am using Windows 10 (64bit).
https://www.giss.nasa.gov/tools/panoply/download/


texts: Joining attributes via webmap hyperlinks

texts:
I have a parcel dataset from Henrico County Virginia but it doesn't include parcel owner names, which would be incredibly helpful in the analysis I'm running. The county has an online parcel webmap and in the attributes it includes a "Real Estate Data" link,

which takes you to a webpage of the county that provides the owner name.

I have the GPIN attribute in my parcel dataset, which could act as a join field/unique ID.
Is there any way I can obtain the owner names and join them to my parcel dataset or to obtain this parcel dataset with the owner info by other means?
My thought is that this would entail some sort of web scraping, which I'm not familiar with. I know basic python but that's as far as my coding experience goes.
I'm using ArcMap, and I often obtain data from state and county Map and Feature Services.


texts: Digitizing from OpenLayers Plugin in QGIS, or create raster?

texts:
Is there an easy method for digitizing from a layer added from the OpenLayers pluggin. As i cannot snap to it or use it as raster so the entire map does not load?

two separate issues there.


texts: Creating ArcSDE Geodatabase without post installation?

texts:
Is there any way to create an ArcSDE SQL Server 2008 geodatabase without :

SDE post-installation execution on the server
having sde service running on the server

and using :

execution of some DDL scripts to create database on the SQLServer
execution of "sdelayer -o register" command to register some
tables (containing Geometry type
fields) as feature class.


texts:
Hello I have created a table with a column type PATH as above (following the Postgres documentation)

CREATE TABLE public.geometry_polyline_volume
(
    id          serial not null primary key,
    distance    float not null,
    height      float not null,
    coordinates path not null
);

Trying to insert the above values as a row

INSERT INTO public.geometry_polyline_volume(id, distance, height, coordinates)
VALUES (2,  500, 0, path((15.878137629895164,47.08306448089695), (15.56169808311181,47.219041634920686), (15.267442604782124,47.4201665137259), (15.092631384557304,47.71366328136526), (15.234428926980286,47.95865145177352)));

I am taking the following error

  Message    : ERROR: column "coordinates" is of type path but
  expression is of type record


texts:
I have the following question about formatting a string which I use to create output file names.

I loop a SearchCursor to create .png file. I extract two columns form the cursor e.g. "GMS_ID" and "ID". (codepart - with arcpy.da.SearchCursor("dbk_object_test_SpatialJoin2",['GMS_ID','ID','SHAPE@']) as cur:)

Then I concatenate the two columns like:

vakkrt_name = str(row[0]) + "_" + str(row[1])

Almost at the end I use:

outpng = r"D:\\arcpy\\output\\DBK_" + format(vakkrt_name) + ".png"

to create the file reference and after that I use:

  arcpy.mapping.ExportToPNG(mxd, outpng,df,
                          df_export_width=1600,
                          df_export_height=1200,
                          resolution=200)

to really create the file on disk.

The result is something like:

DBK_1234.0_5678.0.png

Can you help me by disabling the decimal point probably by using the Format method.


texts: ArcGIS JS, use only right half of the map area when setting extent?

texts:
I have a full page arcgis js map (v3.15). On the left hand side of the map, I have a div panel positioned on top of the map that will contain settings/configurations for a user to choose from (see attached screenshot).

The problem is, when I set the maps extent to match that of a feature layer, a lot of the graphics get covered up by the panel on the left. Is there any way to only use the right half of the map to calculate the zoom and center position so all of the graphics in a layer appear on the right half of the screen?


texts:
Using QGIS v 3.4.4.

I have a polygon layer in which some shapes have a very narrow extension that i would like to remove. Unless zoomed in these extensions appear as dangles but when zoomed right in you can see that they are part of the polygon shape. I have tried the v.clean toolset without success and snap geometries to layer which ends up crashing the project. I can manually edit the shape to correct but I don't know how many or where these occur in the layer.


texts: What are the best automated ways of placing the nodes and edges of a routing graph using the GeoJSON data of a multi-level building? What other data is needed (area layout, accessor directions etc)?


texts: How can I make ArcGIS 10.2 overwrite folders instead of pasting a new/copy folder (with a (1) suffix?)

texts: Just upgraded to 10.2 from 10.0 and I have always edited data on the C: drive then pasted this updated folder onto the "work" server.  In 10.0 the folder is overwritten. In 10.2 the folder is copied and given the (1) suffix giving me two folders.  How can I fix this?  Thanks!


texts: How to implement language in geonetwork?

texts: How can I implement another language in Geonetwork 2.10.4? I have tried to follow this tutorial https://trac.osgeo.org/geonetwork/wiki/HowToTranslate, and in the section where are  described ISO19139 / Multilingual metadata and Translation of metadata schema labels, there is no such path as given in example. I have tried with changing code and haven't managed to fully impelemnt language. Is there any other way?


texts: How to add virtual (e.g. $area) fields to multiple layers

texts:
I wondered if I could add virtual fields to the feature tables of a number of polygon layers in a QGIS 3 project (or a geopackage).
As an example, I would like to add the virtual field area (km²) as round($area/1000000, 2) to ten polygon layers stored in a geopackage.
Is there a QGIS way of doing it? If not, any leads if I do this from R? (I don't use Python.)


texts: Delete all identical records from each set of identical records

texts: I have intersected a point and a polygon layer and I would like to remove all cases where two or more points intersect with the same polygon. I do not want to keep one row of each duplicate, but remove them all, so I cannot use Delete Identical. I have FID:s of the polygons in one column and would like to remove all cases of the same FID occurring two or more times. I am using ArcMap version 10.5.


texts: Ordnance Survey MasterMap Symbols not appearing properly styled on QGIS

texts:
Using plugin OS Translator II on the Exeter sample database.
Ticked the boxes "Add OS Styling Fields" and "Apply Default OS style"

Connected to PostGIS database. All worked apart from the Symbols for CartographicSymbol and TopographicPoint. The symbols are random letters and shapes. And huge. Big Magenta N and Grey question marks for NonConiferenceTrees etc.

I've got the osmmsymbol folder in the appropriate QGIS folder. Changed Settings-Options-SVG Paths to C:/OSGeo4W64/apps/qgis/svg/ in addition to the default C:/OSGEO4~1/apps/qgis/./svg/ (What is that weird dot for??).

I can manually change the symbols. Via LayerProperties-Symbology and selected one of the 6 Legends. A double click to show Symbol Selector and then click on Font Marker and change to SVG marker. But thereafter it is hunt the symbol. And when applied the size is wrong. After fiddling I can fix that too.

I believe the symbols should automagically appear without manual intervention.

Any ideas?


texts: overlaying raster layer hides my vector layers in PyQGIS

texts:
I loaded vector layers and rendered them without a problem but when I tried to add a raster layer as well the map suddenly only shows the raster image and my vector layers are missing.

Anything I can try to diagnose the issue here.  I was following the routine steps in the PyQGIS 2.0 cookbook.

    import sys
    from PyQt4.QtCore import *
    from PyQt4.QtGui import *
    from qgis.core import *
    from qgis.gui import *
    import sys
    import os
    from mapui import Ui_MainWindow

    os.environ['QGISHOME'] = 'C:\\Program Files\\QGIS Brighton\\apps\\qgis'

    class ShapeViewer(QMainWindow, Ui_MainWindow):
    def __init__(self):
        QMainWindow.__init__(self)

        #Required by Qt4 to initialize the UI
        self.setupUi(self)

        #set the title for the app
        self.setWindowTitle("Map Viewer")

        #create the map canvas
        self.canvas = QgsMapCanvas()
        self.canvas.setCanvasColor(Qt.white)
        self.canvas.enableAntiAliasing(True)
        #self.canvas.useImageToRender(False)
        self.canvas.show()

        self.clickTool = QgsMapToolEmitPoint(self.canvas)

        self.layout = QVBoxLayout(self.QgsMapCanvas)
        self.layout.addWidget(self.canvas)

        layers = []

        uri = QgsDataSourceURI()
        #underlying map
        raster_layer = QgsRasterLayer("d:\\working\\hyp_lr\\hyp_lr\\hyp_lr.tif", 'geography layer')

        if not raster_layer.isValid():
            print "raster layer is invalid!"
            return

        QgsMapLayerRegistry.instance().addMapLayer(raster_layer)
        self.canvas.setExtent(raster_layer.extent())
        cl_raster = QgsMapCanvasLayer(raster_layer)

        layers.append(cl_raster)

        uri = QgsDataSourceURI()
        uri.setConnection("my_linux_box","5432","mydb","me","mypwd")
        uri.setDataSource('public','some_table','geom','')
        vector_layer = QgsVectorLayer(uri.uri(), 'Some Data', "postgres")

        if not vector_layer.isValid():
            print "Vector layer is invalid!"
            return

        QgsMapLayerRegistry.instance().addMapLayer(vector_layer)

        self.canvas.setExtent(vector_layer.extent())

        vector_canvas_layer = QgsMapCanvasLayer(vector_layer)
        layers.append(vector_canvas_layer)

        ####
        self.canvas.setLayerSet(layers)

        self.canvas.show()

def main(argv):
    #create Qt application
    app = QApplication(argv)

    #init qgis libs
    qgis_prefix = 'C:\\OSGEO4~1\\apps\\qgis'
    QgsApplication.setPrefixPath(qgis_prefix, True)
    QgsApplication.initQgis()

    #create the main window
    wnd = ShapeViewer()

    #Move the app window to upper left
    wnd.move(100, 100)
    wnd.show()

    #run!
    retval = app.exec_()

    #exit
    QgsApplication.exitQgis()
    sys.exit(retval)

if __name__=='__main__':
    main(sys.argv)


texts: Removing saved custom colors from "favorites" in ArcGIS Pro

texts:
This is easy to do in ArcMap via style manager but I have spent half an hour trying to figure this out to no avail in ArcGIS Pro.
I accidentally saved some duplicates of custom colors and it says the colors themselves were saved in "favorites." Managing that folder only shows a saved polygon symbol and none of the colors I created.
How do I remove the colors I've saved?


texts:
Is this the right way to get centroid x and y of a polygon. How do I get centroid x,y, of a polygon

layerName = "INTERIOR"
layer = QgsMapLayerRegistry.instance().mapLayersByName(layerName)[0]
print layer.name()

for f in layer.getFeatures():
    print f.geometry().centroid().asPoint()


texts: My view in Oracle when displayed in ArcMap/ArcCatalog values are being dropped to 0

texts:
I created a view in an Oracle Database in order to match formatting of a new database/table I am going to be migrating the data to.

The problem is when I look at my view in Oracle I see the correct values (at least values that I expected to see). However, when I bring the view into both ArcMap and ArcCatalog through the add data or database connection those values turn to zeros.

Does anyone have any thoughts?
I am using arc10 and SQL Developer


texts: Calculating column mean of array or list using Google Earth Engine

texts:
I want to calculate the column mean of the following set of data as the result of 1 row and 7 columns.
Which functions should I use in Google Earth Engine or should I write an extra function?
var data=[['1','2','3','4','5','6','7'],
              ['8','9','10','11','12','13','14'],
       ['15','16','17','18','19','20','21']];


texts: Get vector layers of self.iface.legendInterface()

texts:
I have an application which I did not use qgisMapComboboxes; I used simple QComboBox.

without changing the combobox type, how can I get only vector layers in VectorLayers = iface.legendInterface().layers() before my iteration and combobox layer manipulation?


texts: Tool that modifies number of vertices but not the feature itself

texts:
Do you know the reason, why some ArcGIS tools (like Clip and Intersect) modifies number of vertices, even if the feature is not modified itself? Like on the picture, I use Clip and Intersect on line with polygon. The line has 26 vertices before executing function. After, the line has 517 vertices even the tools don't modify that line.

I was only wondering how this is possible and if this is some bug or there is some reason why this tools do it. It happens with both functions (I tried only this two) no matter it's shapefile or geodatabase.

EDIT: The line was created by r.flow algorithm in GRASS GIS, then exported into shapefile. I think it hasn't true curves. The polygon is simply buffer around the line at 10 meters. Coordinate system of data and data frame is EPSG 5514. I'm using Win 7 (64-bit).


texts: Qgis -> DXF2ShapeConverter -> Inscriptions in Qgis

texts:
Sorry, but this is my first question which very often was asked . It has not yet answered very well.
I would like to import a DXF with captions in QGIS. I have already tried this with the DXF converter but it always crashed GIS .
Now my question , are there alternatives where I can get the label of the DXF with in the GIS ?
 Mfg Malte


texts: I am currently trying to figure out how a set of points were plotted without full XY coordinates. I was given Y coordinates, but not X and there are points plotted on a map. I think these points were just placed by drawing points/creating features. Is there any way I can figure out the X coordinates from the Y?


texts: How do you create and use a SQL Server spatial view in ArcMap?

texts:
All,

(using SQL Server 2012, ArcMap 10.1)

I have created a SQL spatial statement to test out how a) the union function works, and b) how creating a view of this data will work with ArcGIS.

a) the function i'm testing is:

SELECT GEOMETRY::UnionAggregate(shape)
FROM CENSUS_TRACTS_2010_2877
WHERE PCT_ASIAN > 5

...which gives the following correct result:

b) our database would greatly benefit using a spatial view of our boundary data (using a union, etc.) to display boundaries to our internal users as a single table.

To accomplish this, it seems that a spatial view would be the best use case, as that gives us the option to edit/maintain one file, and the view will serve up a fresh version of the boundaries.

Problem is, when I create the view in SQL Server, it will not show up in ArcMap, and I get the error:

Warning: Could not add the specified data object to the map.

I think I might be creating the view incorrectly, as this is the field properties ArcMap sees from the view:

...which doesn't seem correct.

Can anyone point me to a place that describes how to create / register / use a spatial view in ArcMap / SQL Server 2012?

Thanks!

-mb


texts: latitude and longitude coordinates into postal code with town names

texts:
Here is my problem.

I have a csv file, with latitude and longitude coordinates. I would need to know
the name of the French towns and the relevant postal codes that match these coordinates.

My initial coordinates seem to be quite precise as I'm using Maxmind Geolite city. I'm trying to match the coordinates with the data set from this website.

http://public.opendatasoft.com/explore/dataset/correspondance-code-insee-code-postal/?tab=export&location=3,18.54017,-3.01253&basemap=mapquest

However, the coordinates when doing the matching does not work properly.

Question is the following: Is there a online tool where I can load a csv file with the coordinates and whave in return, the postal codes and the names of town?

I'm insisting for French towns as I'm in France.


texts: Conversion of raster features to polylines: skeletonization?

texts:
I would like to identify linear features, such as roads and rivers, on raster maps and convert them to a linear vector features (i.e. polylines).  I typically use R as my primary GIS, and have posed this question as it pertains to R, here, at StackOverflow without resolution.  Apologies for cross-posting, however, here I am seeking any means to do this efficiently, (preferably, but not necessarily, in a manner that can be automated using code.)

Converting linear features of a certain class on a raster to polygons can be done using rasterToPolygons() in the R package raster, for example.  But how can I turn these polygons into polylines objects (i.e. skeletonize these polygons to a centre line) efficiently and reproducibly without resorting to to manual steps in vector graphics software, for example.


texts: Join Attributes by Location: .mdb and shapefile in QGIS

texts:
I can join attributes by location in QGIS using 2 shapefiles.

However, I do not know how to do the join when one of my files is in a .mdb format while the other is a shapefile.


texts: Setting transparent to color for raster layer using ArcMap?

texts:
I'm working on a map where I'm overlaying a topographic map on a hillshade map (from ALS data). With the right levels of transparency, those maps tend to be the best base for all kinds of stuff.
However, it appears that my topographic map (that already includes some very light hillshading) has areas that should be white, defined as transparent. So at those, the hillshade data below will always stick out in their pure gray.

(shown here with the topo map with no transparency, to highlight the issue)
How do I set transparent to be white?
Most stuff I find online, is about the reverse, setting a color to transparent, and when it's the other way round, it appears to not work with my data, e.g. the method shown at Setting multiple colors transparent in raster map in ArcMap?
Perhaps it also depends on what kind of raster I have.


texts: Very different static gps tests for different devices

texts:
I don't understand the results that I got when I measured with three different devices that were lying on the ground for couple of hours.
These devices were: Nexus 10 tablet, Samsung Galaxy Core Prime smartphone and Samsung Galaxy Tab (quite old tablet). I can add more info if needed.

These devices were lying right next to each other (which means that the receivers were few centimeters from each other). I used the app GPS Logger (it has red symbol in Google Play) to capture the coordinates for couple of hours.

Here is the graph with all the points:

So my question is, why did Nexus 10 (the middle blue color) measure with such an error shooting points everywhere and constantly shifting while Core prime (light blue) was moving slowly over couple of hours just wandering around and Galaxy Tab (dark blue) was just switching points after tens of minutes with no path whatsoever? Why is there such a difference and what causes it?

Other info: it was a clear day and I measured for over 6 hours.


texts: Sentinel 2 Atmospheric Correction Inconsistencies

texts:
I'm having some peculiar issues with Sentinel 2 data. I've downloaded two coterminous tiles (one L1C tile from EarthExplorer and one L1C/L2A tile from AWS) taken at the same time and date. I have done this to verify that my methods of atmospheric correction (ESA's sen2cor tool) are valid and give me precisely what I would get from simply downloading AWS' L2A product. The tile ID follows: S2B_MSIL1C_20190113T181719_N0207_R084_T12TWR_20190113T200650

The issue is, it doesn't produce identical results. When I go to calculate NDVI from my bottom-of-atmosphere product, I get different values when using AWS' L2A versus the L2A I produced using sen2cor.

I've done a little research on this and found reference to a "pilot" introduced to ESA's L2A_GIPP process in 2018. There seems to be a version 2.11 floating around out there, while the process in sen2cor (2.07) produces "L2Ap" instead of "L2A"

My question is: Does anyone know why an identical atmospheric processing chain would produce dissimilar results?


texts:
I am using FTTx Planner (Open Source element of Geospatial Network
Inventory system) but when I try to use GniPlannerFTTH plugin in QGIS I get the following error:

  "Couldn't load plugin GniPlannerFTTH due to an error when calling its
  classFactory() method"


texts:
I just want to draw rectangle for style.
if I do this

 image: new ol.style.RegularShape({
        fill: fill,
        //stroke: stroke,
        // offsetY: -25,
        points: 4,
        radius1: 30,
        //radius2:20,
        angle: Math.PI / 4
    })

it gives me this square

else if i do this

 image: new ol.style.RegularShape({
        fill: fill,
        //stroke: stroke,
        // offsetY: -25,
        points: 4,
        radius1: 30,
        radius2:10,
        angle: Math.PI / 4
    })

it gives me this

How can I draw a rectangle?


texts: How can I make a buffer with xxx meter by NetTopologySuite

texts: I am very new in NetTopologySuite, I work with the coordinatesystem in WGS84 by lon/lat, with the NetTopologySuite. I found the Buffer function to create a buffer with a distance, but how can I make a buffer with distance by the unit of meter? e.g make a 200 meter buffer of one point.


texts:
I have two tables, A and B.

I want to Union(/Dissolve) B and then difference A,B. To get whatever is left of A when overlapping B is erased:

select a.someid, st_difference(a.wkb_geometry,st_union(b.wkb_geometry))
from a, b
where b.group_id=43 and b.type_id=12 and st_intersects(a.wkb_geometry, b.wkb_geometry)

I get an error on line 1:

  ERROR: column "A.wkb_geometry" must appear in the GROUP BY clause or
  be used in an aggregate function

How to fix this? Without st_union I sometimes get duplicates a's (on id, not geometry) in the output, which I don't want. Maybe it is better to Union the output?


texts: Kriging for unevenly spaced points using R?

texts:
I started to read about Kriging and the implementation into R. In all the tutorials I found, a (regular) grid is used to do ordinary or universal kriging. My situation is as follows, I try to interpolate house prices, which of course do not lie on a regular grid. For all houses I have a set of covariates including features such as number of rooms and space. For a subset (approx. 20%), I have the corresponding house price. When doing Regression-kriging in estimating price ~ rooms + space, I cannot provide the data for the regular grid points.
Is there a way of using unevenly spaced prediction points or how is this best tackled?


texts: Visualize text file with location info and intensity

texts:
I have an ascii text file containing location data (column 9-lat and 10-long) and intensity(column 20)

  200010  207  020311    40658.5  406593 52 344927.31 7100203.50  -26.2078720  127.4491855 345060.64 7100369.14   26.4  650.3  628.0 55471.293    20.168 55648.817 55637.523  -146.062

the text file has many lines 10k+

I am trying to export this ascii text into a GeoTIFF using GDAL, but not sure how to proceed.

Ideas?


texts:
Since I am a beginner to postgis I am facing problem with finding points within polygons from my table.

This is my table structure

   Column    |          Type           |
-------------+-------------------------+
 id          | bigint                  |
 type        | bigint                  |
 geoproperty | geometry(Geometry,4326) |

types are 1=point,2=road,3=polygon

I want to get all the points within a specified polygon.

Can somebody help me to get a working query?
I tried using St_within and st_contains but failed to get results.

select id from test_table g where id = 123 and ST_Within(g.geoproperty,g.geoproperty);


texts: NED (LTP) to ECEF Orientation Conversion (Tait Bryan Angles)

texts:
I have an application where I need to convert an orientation from NED - given a WGS84 latitude and longitude - to ECEF. I have yaw, pitch and roll as Tait Bryan angles (XYZ) relative to the local tangent plane but I need this in geocentric format.

Initially I was reading a few white papers on ECEF coordinate conversion but all those I could find only seem to cover position and heading. I'm able to calculate both position and heading correctly but I'm not sure what is is involved to transform pitch and roll. Roll is particularly tricky as a standard 3D velocity vector doesn't include a roll component.

My gut feeling is the solution is going to involve converting to quaternions or DCM and coming up with an appropriate transformation matrix but I haven't been able to find any examples of this being done. Any advice would be appreciated.


texts: How can I select ways that have missing access parameters "[key]=Not Specified"?

texts:
Follow up to How can I programmatically identify missing sidewalks in OpenStreetMap data?. I've noticed that a number of streets and roads in a region I'm interested in are missing pedestrian access when they have sidewalks, and this is affecting pedestrian Mapbox directions.

How can I find ways in Overpass Turbo that are missing this access parameter?


texts: Way to get an adjusted attribute table after clipping raster with polygon?

texts: I'm relatively new to ArcGIS in general, so please excuse any incorrect terminology. Currently, I am trying to clip raster data using polygons. I have clipped six of these, but the attribute table that resulted from five of them is just the attribute table for the original raster data. Strangely enough, the sixth one has adjusted values and is entirely usable. The solution I've seen for this is to re-build the attribute table, but this gets rid of attributes that I've manually added. Would anyone be able to offer advice on either what happened for that one that is correct, or otherwise on how to adjust the table without deleting my manually added attributes?


texts: Guidelines for publishing large datasets on QGIS Web Client

texts:
Can anybody provide simple guidance to publishing large datasets on QGIS web client? I have a mixture of shapefiles and rasters to display. Biggest shapefile is  230 MB. Biggest raster is 1.4GB.

At present the browser shows a loading icon that does not go away...

Do I need to be creating pyramids or is this handled by the server? Maybe there are other quick solutions people can suggest....


texts: There are several ways to calculate slope in qgis using a DEM layer, but how does it actually work? My results show pixels with a slope of 40°+, this seems like a lot (in my DEM, one pixel is 25*25 m). Does it mean that the average slope of the pixel is for example 30° and in reality it could be that the pixel is actually quite flat, but there's a small part of it with really steep slope?


texts: An expression to selection within subsets of vector layer in QGIS

texts:
Does anyone have an idea on "an expression to selection within subsets of vector layer in qgis"? Qgis has “Random selection within subset” to select randomly but I have to select within subset with the expression where a condition is:

WHEN case
there are more than one class, then select only one '1',
if there is no '1' then select '2',
if there is no '1' & '2' then select a '3'.

an


texts:
I think I've run into a bug in how SLD rules are applied in geoserver. This is on GS 2.12.2 with the css styling extension. In my table I have units that I'd like to render with different icons and different colors based on a few attributes. I'm using red-blue-green cars-trucks-planes. I'd also like to show labels depending on zoom level. I'm testing this with cars only. Currently, when I'm zoomed out only red cars are being rendered, with no labels. The green and blue icons aren't being rendered until I zoom in past 20M, at which point the labels also show up. If I remove the '[@sd < 20M]' from the label rule, everything is shown as expected. I don't entirely understand the ordering in which sld rules are applied, but the fact that only some mark icons are being shown makes this feel like a bug. Let me know if you need any more info.

@mode "Flat";

/* @title Truck Red */
[status = 'A' OR status = 'B'] [vehicle_unit_detail_id IS NOT NULL]  {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Truck Green */
[status = 'C' OR status = 'D' OR status = 'E' OR status = 'F'] [vehicle_unit_detail_id IS NOT NULL]  {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Truck Blue */
[status <> 'A'] [status <> 'B'] [status <> 'C'] [status <> 'D'] [status <> 'E'] [status <> 'F'] [vehicle_unit_detail_id IS NOT NULL]  {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Car Red */
[status = 'A' OR status = 'B'] [vehicle_unit_detail_id IS NULL] [service = 'P'] {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Car Green */
[status = 'C' OR status = 'D' OR status = 'E' OR status = 'F'] [vehicle_unit_detail_id IS NULL] [service = 'P'] {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Car Blue */
[status <> 'A'] [status <> 'B'] [status <> 'C'] [status <> 'D'] [status <> 'E'] [status <> 'F'] [vehicle_unit_detail_id IS NULL] [service = 'P'] {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Plane Red */
[status = 'A' OR status = 'B'] [service = 'N'] {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Plane Green */
[status = 'C' OR status = 'D' OR status = 'E' OR status = 'F'] [service = 'N'] {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

/* @title Plane Blue */
[status <> 'A'] [status <> 'B'] [status <> 'C'] [status <> 'D'] [status <> 'E'] [status <> 'F'] [service = 'N'] {
  mark: url('http://localhost:8080/geoserver/www/styles/images/...');
}

[@sd < 20M] [@labels(true) = true] {
  label: [vehicle_id], [location_name];
  label-anchor: 1.5 0.5, 1.2 2;
  font-family: "Trebuchet MS Bold";
  font-fill: white;
  font-style: normal;
  font-size: 12;
  halo-radius: 0.5;
  halo-color: black;
}

EDIT1: I've tried moving the labeling rule into its own FeatureTypeStyle, but I got the same behavior.

EDIT2: I did a small test to see if the issue had something to do with location, but when I created 2 objects in the same location, the red car rendered when zoomed out and the blue one rendered over it when I zoomed in. Doesn't seem like a tiling/location issue.

EDIT3: Here's a pastebin of the generated sld if anyone wants to take a look. https://pastebin.com/cRgDDFGs


texts: GeoServer and PostGIS connection technology?

texts: Does anybody know how GeoServer connects to PostGIS. Which technology is used? Is it JDBC?


texts: Publishing PostGIS raster vs publish file raster in GeoServer; what's best regarding performance?

texts:
I want to show high resolution raster layers in a web app (Open Layers) using GeoServer and I need to know what option is best in terms of performance.
My ideas are:

Publishing the raster layer as a file in GeoServer
Add the raster in PostGIS and use GeoServer to extract the raster from the db table.

What are the advantages and disadvantages of the options, or which other options should I consider?


texts:
I am trying to make a QGIS plugin to change which band in a multi-band data set (netCDF where the band represents time steps) is used for a gray scale rendering. The file works fine with manual changes, but I cannot find how to change it programatically. I have set up a dialog box where I can select which layer to work on:

for layer in self.iface.legendInterface().layers():
        # iface.messageBar().pushMessage("Type",layer.name(), level=1, duration=10)
        if type(layer) is QgsRasterLayer and layer.bandCount()>1:
            self.layerlist.append(layer)
            self.dlg.comboBoxLayer.addItem(layer.name())

my QGIS related imports are (as far as I can see)

from qgis.core import *
from qgis.utils import *

I have found the attribute setGrayBandName that should select which band to use for presentation, but when I run (knowing 0<=sldValue <=number of bands)

self.dlg.labelStartTime.setText(str(sldValue))
self.activeLayer.setGrayBandName(self.activeLayer.bandName(sldValue))
self.activeLayer.triggerRepaint()

it crashes with "AttributeError: 'QgsRasterLayer' object has no attribute 'setGrayBandName'"   - if setGrayBandName does not exist, what should I use then? And why do I find a number of references to it, but none saying that is has been taken out of use? - or have I just misunderstood something? Do I need to import som more modules?

seems like I misunderstood something, I need to set a new gray band on the renderer, i.e.

self.activeLayer.renderer.setGrayBand(sldValue)

But, this only works if I start in the console by selecting the layer and then doing

layer=qgis.utils.iface.activeLayer()
rndr=layer.renderer()
layer.renderer=rndr

Then I can start playing around. If I try to do the rndr=layer.renderer() in the code in the plugin, I am told something like that type of object is not callable...

The present version of the module is available at https://github.com/sickel/TimeStepSelector


texts:
I have coordinate system with x and y. Also I need to show there x,y point but I have coordinates latitude, longitude.

How can I transform latitude, longitude into x,y to show on x,y scale?


texts: "Spatially normalizing" rasters of different extents to sum values of overlapping cells in ArcGIS Desktop?

texts:
How do I add rasters together that are of different spatial extents?

I have many, many features that overlap over a large area (entire east coast) and I am looking to sum them "stacked" where they overlap. My strategy is to:

1) Convert each feature to a raster (using Polygon to Raster tool)

2) Set Null values for feature to 0 using Con(IsNull("raster"),0, "raster") in Raster Calculator

3) "Normalize" each raster (spatially) to have matching spatial extent (of the entire study area) - do this by adding each raster to the "Constant Raster" I created, which has the study area's extent, the same projection, and the same cell size of the individual rasters. The result should be a new raster to the extent of the study area, with the same correct values at the location of the converted-feature in question.

4) Use Raster Calculator to simply add the normalized rasters together.

Something's not working at step #3. When I add the projected, constant raster (all values = 0) to the individual feature-converted-to-raster (with null values set to 0) in the Raster Calculator Tool, with Processing Extent set to that of the Constant Raster, the resulting raster output is of the correct extent but the values of the feature-converted-to-raster are incorrect. Instead of an expected range of, for example, [0, 77] (where the feature that was converted had the value 77, and any "background" value = 0) the value is something extreme, such as [-1.79e+308, 1.79e+308] [Edited to explain unexpected output]

Can someone explain what it is I am missing here? My searches on StackExchange only revealed this question but I must have missed something.


texts: Getting Landsat NDVI data processed with LEDAPS on USGS earthexplorer

texts:
I am trying to download NDVI values from Landsat on the earthexplorer. I read there were higher-quality NDVI series processed with the so-called LEDAPS algorithm. The page from the LEDAPS website mentions:

  To obtain Landsat-5,7 surface reflectance products processed with
  LEDAPS atmospheric correction, you may now order them directly from
  USGS via Earth Explorer: (select the "Landsat CDR" tab under Data
  Sets).

I do not see however any Landsat CDR tab on earthexplorer?

How do I download NDVI data from Landsat on earthexplorer?

Is that dataset atmospheric-corrected?


texts: How to update a line attribute in PostGIS using an attribute from the polygon in which it is contained?

texts:
I have a road dataset which I want to update using the speed information from the state polygon in which it's contained.

for all roads (line) within a state (polygon)
if the road speed is null
road speed becomes state speed

My preference is to use Contains or Within and only attribute speed information to the road if the entire line geometry is completely inside a polygon.

How do I do this in PostGIS?


texts:
my idea is to add different buttons for drawing interaction in OL3 such as :

<button id="Point" onclick="addInteractions(Point)">button1</button>
<button id="Polygon" onclick="drawPolygon(Polygon)">button2</button
<button id="LineString" onclick="drawLine(LineString)">button2</button>
var drawstyle = new ol.style.Style({
          fill: new ol.style.Fill({color: 'rgba(255, 255, 255, 0.2)'}),
          stroke: new ol.style.Stroke({color: '#ffcc33', width: 2}),
          image: new ol.style.Circle({radius: 7, fill: new ol.style.Fill({        color: '#ffcc33'}) })
      });
var source = new ol.source.Vector();
var vector = new ol.layer.Vector({source: source, style: drawstyle)};
function addInteractions(str) {
        var draw = new ol.interaction.Draw({
          source: source,
          type: str
        });
        map.addInteraction(draw);
}

on the official site of OL, all examples used the html select tag, but I don't want to use it.
can someone help me please !


texts: Finding pixel or map size of TextSymbol using ArcGIS Runtime SDK for .NET?

texts:
I'm using the ArcGIS runtime SDK for .Net, version 10.2.7.

I have a point graphic that is symbolised with a TextSymbol.  I would like to add a separate polygon graphic to the graphics overlay that covers the extent of the text symbol.

What I have tried is below.  It works fairly well, but there is a small difference between the text and the polygon.

 var symbol = GetGraphicSymbol();
 if (symbol == null)
     return;
 var text = symbol.Font;

 var size = graphics.MeasureString(symbol.Text,
     new Font(text.FontFamily,
              (float)text.FontSize,
              FontStyle.Regular,
              GraphicsUnit.Document));
 if (size.Width == 0 || size.Height == 0)
     return;

 var width = mapView.UnitsPerPixel * size.Width;
 var height = mapView.UnitsPerPixel * size.Height;

I understand that on any scale change the polygon graphic would need to be recalculated.


texts: Counting records that match attribute from another feature class?

texts: I have geometry feature class of states that I would like to build a choropleth map from by using the quantity of people in that state. However, the data of what people in what state is in seperate table with other attributes like first and last name and so on. Is there a way to count all of the, say, Washington State people from the seperate table and append a "population" attribute at the end of the geometry feature class' attribute table so that I may then use this newly derived "population" attribute to use for the choropleth map?


texts: Calculating the percent of binary values in grid of 1 km by 1 km across the entire raster

texts:
I have a raster of the entire united states with 30 meter by 30 meter pixel size.  Each pixel has a value of either 0 or 1, with 1 defining a certain type of landcover.
I need to find the percent of 1's within 1 km by 1 km areas gridded across the entire raster so that I know for each 1 sq km, what percent of it is the landcover I'm interested in.
Despite a lot of googling, I am unsure about how to proceed.  Mostly how to create the grid.  It looks like once I have that I could use v.rast.stats in GRASS.  But maybe I'm way off.
As a result, I don't have a specific error I'm encountering, but rather do not understand the tools well enough to know how to proceed to get my desired output.


texts: Clickable HTML link in QGIS print composer pdf export?

texts:
According to qgis help files here:

  'Labels can be interpreted as html code: check the checkbox Render as HTML. You can now insert a url, an clickable image that link to a web page or something more complex.'

And yet, my label with the contents:

google

produces a pdf that's rendered with 'google' in blue text like a link but it's not clickable.

Is it possible to have clickable links in the output? And if so, what am I doing wrong?


texts: How to add style to GeoServer through REST API

texts:
Let imagine. the style content is like below CSS code.

[status = 'Layih&#x259;']{
  fill: #e1e1e1;
}

[type = 'Sahil'] {
    fill: #6f9bcc;
  [@sd < 6k]{
    stroke: #3235ff;
    fill: url("https://imgur.com/BdQQhxA.png");
    fill-mime: "image/png";
  };

}`

I can not add the content of style to geoserver. It only creates a style name and style file name. Could you please help?

   await axios({
        method: 'POST',
        headers: {
            'Content-type': 'application/json',
            Authorization: 'Basic YWRtaW46Qm94eUFzdHJpZDYzMTk=',
        },
        url: 'http://localhost:8080/geoserver/rest/workspaces/sde/styles',
        data: {
            style: {
                name: 'style_clone',
                filename: 'style_clone.css',
            }
        }

    })
        .then(
            response => {
                return res.status(201).json(response.data);
            },
        ).catch(error => {
            console.log(error)
            return res.status(500).json(error);
        });;


texts: Automating calculation of Population Density in Census Tracts?

texts:
I am quite new to Python and I am attempting to practice by making a simple program with ArcGIS. I wanted to make a script that automated the calculation of population density of selected census tracts within the city of Plano, TX. I also wanted to give the ability for the user to choose the census tract.

I then wanted to provide a map layout in ArcMap of the specific tract with the population, area, and population density being shown.

I realize most of this can be done using ModelBuilder, but I needed practice with Python and thought this would be a good simple program to start with. Anyone have any pointers on the script and also how I could possibly make it better?


texts:
I ran las dataset, then las to raster using max for DSM.  Two and a half hours later I have 3 areas showing pixilation using higher level coloring.  See enclosed image.

I found on web that seamless might work though 2,2 might be small.


texts:
I have a batch of (binary) ESRI files containing 25m resolution raster data.  I need to check my interpretation of the OSGB coordinates therein.

gdalinfo reports the following for a file containing a 3578 x 3833 raster:

Pixel Size = (25.000000000000000,-25.000000000000000)
Corner Coordinates:
Upper Left  (  486425.000,  222650.000) (  0d44'37.94"W, 51d53'42.91"N)
Lower Left  (  486425.000,  126825.000) (  0d46' 2.31"W, 51d 2' 1.48"N)
Upper Right (  575875.000,  222650.000) (  0d33'18.65"E, 51d52'27.15"N)
Lower Right (  575875.000,  126825.000) (  0d30'27.19"E, 51d 0'48.00"N)
Center      (  531150.000,  174737.500) (  0d 6'43.71"W, 51d27'21.30"N)

The clarification I need (given the negative Y pixel size) is whether the upper left coordinate given corresponds to the NW corner of the top-left-most cell, or the SW corner.


texts: Python script for creating geodatabase in ArcMap s "Home" directory

texts:
How can I write a Python script that I could add to ArcToolbox which will when run will automatically create a geodatabase in the Home (current MXD file's working directory) and set it as the scratch workspace?

I am new to Python (only real coding I do is VBA) so I've just been lurking and googling for snippets of code to come up with my desired result.  So far I've been stuck on accessing the Home directory; the best I could come up with is a script that will print a message which includes the current mxd filename w/the filepath.  I'd need to remove the filename from the output string to even begin attempting to create a new geodatabase at the desired location.  Below is my attempt at Python coding:

import arcpy

output=arcpy.mapping.MapDocument('CURRENT').filePath

arcpy.AddMessage("The current file is {}.".format(output))


texts:
I'd like to upgrade QGIS from 2.8.1 to the lastest LTR version.

I'm looking for a clear explanation on how to achieve that on Windows (I'm on Windows 7). Such information is not easy to find, including on the QGIS official website that only gives details about installation.

Is it just necessary to download/install the new version or do I have to uninstall programs before that? I don't want to keep the older version alongside the new one.


texts: How to reconstruct a network with multipart edges to a network with singlepart edges

texts: I have a network shapefile that between its vertices , there are edges that have multiple geometric objects. When I try to node it in postgis, it creates multiple segments in the database that topologically are actually one single edge.To clarify, imagine from node A to node B, I have one edge. In the database, this topologically one edge, has been divided into multiple segments and saved in the database. Could someone please let me know how I could reconstruct the network to have a network with singlepart edges in ArcGIS-Desktop or using any other tool?


texts: Problems with saving changes in layerproperties (Qgis)

texts:
I have problems with layerproperties in Qgis. I load a shapefile, edit textsize and color and do labeling some points with text. After i saved the layer as shapefile again and open it......all changes as edited color and textlabeling are gone. Do someone know what exactly the problem could be?

Regards Prospekteur


texts: QGIS - create linking lines or arrows to points when I move lables manually

texts: When I move labels manually in QGIS by using the data-defined placements X and Y fields, I wonder if a line or arrow can be inserted automatically which will link the new label position to the point in question?


texts:
I have 3 geometry tables a, b and c with geometry column geom

What I want is to first check if geom of b and c, intersects under some given condition using a subquery, and feed the subquery geometry resultant rows into another ST_Intersection query, to check whether it intersects the geom of table a as show in the query below.

select geom from a
where st_intersects(a.geom, (
  select b.geom from b, c
  where st_intersects(b.geom, c.geom) = 1
     and b.col_1 >= 180
     and b.col_1 <= 350
     and c.col_1 >= 1800
     and c.col_1 <= 2800)
) = 1

I'm currently getting only result from intersection of the first row of a.geom with the the subquery, how can I get all geometry from table a that intersects result from the subquery of intersection of table b and c under specified condition?