texts: Creating Polyline shapefile from text file using ArcPy?

texts:
I'm trying to create a polyline shapefile from a text file which provides Name, coordX, coordY.  I've created the following script, which I feel will work. However, I can't seem to get it to read the text file properly. I am getting the following error, RuntimeError: Object: CreateObject cannot create geometry from inputs.

import arcpy

arcpy.env.workspace = "C:\\Data"

outFolder = "C:\\Data"
fc = "Paths.shp"

spatRef = arcpy.SpatialReference(26913)
arcpy.CreateFeatureclass_management(outFolder, fc, "POLYLINE", "", "", "", spatRef)

coordinateList = open("C:\\Data\\Paths.txt")

for line in coordinateList.readlines():
    print line

pointList = arcpy.Array()

for x, y in coordinateList:
    point = arcpy.Point(x,y)
    pointList.add(point)

polyline = arcpy.Polyline(pointList)

cursor = arcpy.da.InsertCursor(fc, "SHAPE@")
cursor.insertRow([polyline])

del cursor


texts: Extracting values from multiple raster dataset at location point in one single dbf using ArcGIS Desktop

texts:
Suppose I have 250 raster and my only concern is just a single point. The raster is arranged by ascending date. So if I extract a point value based on the order into one single file then I can graph it just easily.

How do I do it using ArcGIS Desktop?


texts: Publishing Image Services in ArcGIS Server vs. Tile Cache Data Store

texts:
Historically I've published images services through ArcGIS Server, but I now have an ArcGIS Enterprise tile cache data store on a dedicated server and I'm considering moving some image services over to it. I have successfully published a large image service to the tile cache data store and it works great. Our image services are basic cached image services which do not use any mosaic dataset functions, and they are static image sets that do not change. Moving image services to would relieve some of the pressure on ArcGIS Server. So, my questions are these:

Is there a performance difference between ArcGIS Server image services and the tile cache data store?
Are there any negatives to moving image services to the data store?


texts:
What is a ballpark transform?

I read QGIS 3.12 - Setting Project default CRS causes "Ballpark Transform"
and
QGIS warning ballpark transform from
but found the answers not really helpful.


texts: Grass 6 tools missing in Processing tool box

texts: I have installed QGIS via osgeo4w-setup (64Bit). I cannot find grass 6 tools in the toolbox. I also have the dev version of QGIS on my notebook and there I can see and work with the grass 6 tools. Do I have to install something else?


texts:
My Environment is PostGIS, Geoserver, and Openlayers3.

I've a PostGIS table which is published as a Geoserver layer. I need to filter its displayed features according to some complex filters (some complex queries on multiple DB tables). The user would select some criteria from a GUI, which will be used to filter that layer.

What is the best approach to do so?

I've the following (not complete) ideas:

1- Implement the filters as Geotools functions, then call these functions from the Style Document of the layer. But how to pass the user's criterias to these functions?

2- Implement the filters as Geotools functions, then can I call these functions from WFS filters? If yes, there is another problem which is the large size of layer's features, I prefer to display the layer as WMS rather than WFS.

3- Any other approaches?


texts: How to create a projection from a CRS string using pyproj?

texts:
I'm trying to convert between two coordinate systems using pyproj:

import pyproj

p1 = pyproj.Proj("+init=EPSG:4326")
p2 = pyproj.Proj("""PROJCS["NAD27 / California zone II",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982138982]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",39.83333333333334],PARAMETER["standard_parallel_2",38.33333333333334],PARAMETER["latitude_of_origin",37.66666666666666],PARAMETER["central_meridian",-122],PARAMETER["false_easting",2000000],PARAMETER["false_northing",0],UNIT["Foot_US",0.30480060960121924]]""")

p1 creates fine but p2 gives projection not found

I have corresponding points in both coordinate systems but transforming them using EPSG:26742 doesn't give the correct result:

import pyproj

p1 = pyproj.Proj("+init=EPSG:4326")
p2 = pyproj.Proj("+init=EPSG:26742")

a = (-121.330272, 38.547287, 5.661467)
b = (2191840.924977, 321431.770306, 44.012074)

print "Source point", a
print "Target point", b
print 'Converted: ', pyproj.transform(p1, p2, *a) # different :()


texts:
This is my first project in QGIS and I'm also quite new to Python.  PyQGIS is more than a challenge for me!
I found code in How to programmatically set a rule-based renderer in QGIS? to set a rule based renderer which is exactly what I need for my project! I Already slightly modified it (only the rule expressions), so that it would fit into my project.
I got all my "sub points" listed with the correct color in my layer window. But the actual points on the canvas are gone!
Here is my function so far:
def change_color2():
    try:
        inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('points')[0]
        print ('points-layer found')
        idx = inputlayer.fieldNameIndex('amenity')
        if not idx:
            print ('"amenity" does not exist')
        else:
            print ('"amenity" does exist')

            # define some rules: label, expression, color name, (min scale, max scale)
            poi_rules = (
            ('restaurant', '"type" LIKE \'restaurant\'', 'red', None),
            ('school', '"type" LIKE \'school\'', 'green', None),
            ('pub', '"type" LIKE \'pub\'', 'orange', (10.0, 10000.0,)),)
            # create a new rule-based renderer
            symbol = QgsSymbolV2.defaultSymbol(inputlayer.geometryType())
            renderer = QgsRuleBasedRendererV2(symbol)
            # get the "root" rule
            root_rule = renderer.rootRule()
            for label, expression, color_name, scale in poi_rules:
                # create a clone (i.e. a copy) of the default rule
                rule = root_rule.children()[0].clone()
                # set the label, expression and color
                rule.setLabel(label)
                rule.setFilterExpression(expression)
                rule.symbol().setColor(QColor(color_name))
                # set the scale limits if they have been specified
                if scale is not None:
                    rule.setScaleMinDenom(scale[0])
                    rule.setScaleMaxDenom(scale[1])
                # append the rule to the list of rules
                root_rule.appendChild(rule)
            # delete the default rule
            root_rule.removeChildAt(0)
            # apply the renderer to the layer
            inputlayer.setRendererV2(renderer)

    except IndexError:
        print ('points does not exist')

...do I have do take all possibilities of "amenity" into the renderer rules? Because I just picked those of interest by now...
And isn't it necessary to iterate through the features to update them and also important the refresh map canvas? If so, I unfortunately don't have a clue how/where to implement in the existing code.
Some screenshots:


texts:
Installing the Boundless SDK on ubuntu folowing this documentation and getting the below error after running the command sudo apt-get install opengeo-webapp-sdk

E: Unable to locate package opengeo-webapp-sdk

I have not installed any other boundless package. Has the package been removed or is it temporarily unavailable?


texts: How to concatenate values from postgis table which have the same id into a new column?

texts:
I have the following table in postgres-9.3 postgis 2.1. (geom values is too long , I cut to make it short in below table):

id | name | geom
-------------------------
12 | A    | "01030000000"
12 | B    | "01030000000"
12 | C    | "01030000000"
12 | D    | "01030000000"
34 | E    | "050000005B41"
34 | F    | "050000005B41"
34 | G    | "050000005B41"

I would like to group and concatenate all the values into one column where they have the same id number. Expected result as below:

id | all_name    | geom
----------------------------
12 | A,B,C,D | "01030000000"
34 | E,F,G   | "050000005B41"

Could someone please help me?


texts: How to work with Microstation DGN v8i files in QGIS?

texts:
I have a set of multiple files in DGN format (Microstation V8i). I would like to work with them in Qgis.

Since QGIS only support pre-V8i dgn files, I cannot directly open the files.
My first idea was to find a way to publish DGN files on a WMS service and then import it in QGIS.

Do you know any software that will allow me to publish dgn files in WMS ?
Is there a better way to work with dgn files in QGIS?

A little precision : it is not a one time work that I want to do : people work all the times on these dgn files and I need to open them in QGIS.


texts: How to install GDAL and QGIS on a Raspberry Pi?

texts:
I would like to run GDAL QGIS on Raspberry Pi to act as an affordable GIS processing engine.
The plan is to buy about 50 Pis.

Does anybody have any links or information on installing and running GDAL on Raspberry Pi or just on installing GDAL in RASPBIAN environment?

As best I can see Raspbian is my top choice and it is Debian based but ARCH, RASPBMC, PIDORA, RISC OS, and OPENELEC are the flavors of linux open to me out of the box.

Hopefully with Python.


texts: .Length() method from ogr is crashing python

texts:
I created a simple 1-feature polyline shapefile in QGIS 2.18. I would like to now read this shapefile in gdal using python and extract its length.

I have done the following

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
data_source = driver.Open('simple_line.shp')
layer = data_source.GetLayer()
line = layer.GetFeature(0).geometry()
line.Length()

The last line crashes python:

  python.exe has stopped working

  A problem caused the program to stop working correctly. Please close
  the program.

I guess it is to do with my shapefile since when I create a geometry using more "conventional" gdal means I get the following:

import ogr

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(1116651.439379124, 637392.6969887456)
line.AddPoint(1188804.0108498496, 652655.7409537067)
line.AddPoint(1226730.3625203592, 634155.0816022386)
line.AddPoint(1281307.30760719, 636467.6640211721)
line.Length()

The output to the above is 170573.31786365344 which implies it works.

What could be the reason that my shapefile-derived geometry object is crashing python?


texts: Add a PostGIS layer in a QGIS project as a table with no geometry

texts:
As a example, if I check the option Also list tables with no geometry on the Add SpatiaLite layer(s) dialog, it's possible to add to the project a table with a geometry column as a table with no geometry:

Is it possible to do the same thing with PostGIS? (I'm using QGIS 2.18.8)


texts: Accurate area calculation for countries worldwide individually

texts:
I have polygons of countries with its administrative divisions, I would calculate the areas of each division in a country accurately.

I got the bounding boxes for each country, then I divided its width (longitudely) by its height (latitudely).

If the

ratio = 1, the country extends equally South-North as West-East
ratio < 1, it extends West-East
ratio > 1, it extends South-North

But I made my own criteria:

If the ratio < 0.7, I used an Albers equal-area projection +proj=aea with the 1st SP parallel is 1/6 degree below the northernmost latitude of the country, and the 2nd SP is 1/6 degree above the southernmost latitude.
If the ratio is between 0.7 and 1.3, I used Lambert azimuthal equal area +proj=laea centered at the center of the bounding box.
If ratio > 1.3, I used transverse cylindrical equal-area +pro=tcea with central meridian is at the center meridian of the country and k_0=0.9996

My question
Would someone approve it to me, that it lets me get accurate areas relatively better than a projection for the whole continent or the whole world, or isa  it boring time taking method?


texts: Write projected array to netcdf file, best practice

texts:
I'm trying to write projected rasters to netcdf files including projection information in a way that gdal can understand it. I found two ways to achieve that:

By adding the geotransform as an attribute of the netcdf variable containing crs information
By adding x and y (or lon, lat) coordinate arrays as netcdf variables

gdal seems happy with both options, but I'm wondering whether there is any preferred way/best practice? Which option best comply with netcdf conventions?

See below for reproducible examples of both options.

import netCDF4
import numpy as np
from rasterio.crs import CRS
from affine import Affine

file1 = '/path/to/test_0.nc'
file2 = '/path/to/test_1.nc'

# rasterio like raster meta information
meta = {'width': 3000,
        'height': 2000,
        'crs': CRS({u'lon_0': 0, u'ellps': u'WGS84', u'y_0': 0, u'no_defs': True,
                    u'proj': u'laea', u'x_0': 0, u'units': u'm', u'lat_0': 0}),
        'transform': Affine(1000, 0, 0, 0, -1000, 0),
        'count': 100,
}

# Array of random values
array = np.random.rand(meta['height'], meta['width'])

# Util function to generate x and y dimension arrays (only used in second case)
def meta2dim(meta, side):
    """Generate a coordinates array from raster file meta dictionary

    Args:
        meta (dict): dictionary of raster metadata (as returned by rasterio).
            Important keys are height, width and transform.
        side (str): 'width'(lat or y array) or 'height' (lon or x array)

    Return:
        A list corresponding to center row (or columns) coordinates.
    """
    aff = meta['transform']
    if side == 'width':
        res = aff[0]
        start = aff[2]
    elif side == 'height':
        res = aff[4]
        start = aff[5]
    steps = meta[side]
    array = [start + x * res + (res / 2) for x in range(steps)]
    return array

# First, using only geotransform
with netCDF4.Dataset(file1, mode='w') as src:
    # Create spatial dimensions
    x_dim = src.createDimension('x', meta['width'])
    y_dim = src.createDimension('y', meta['height'])
    # Create temporal dimension
    t_dim = src.createDimension('time', None)
    # chlor variable
    chlor = src.createVariable('chlor_a', np.float32, ('time','y','x'),
                               zlib=True)
    chlor.grid_mapping = 'laea' # This corresponds to the name of the variable where crs info is stored
    chlor[0,:,:] = array
    # laea variable to store projection information
    laea = src.createVariable('laea', 'c')
    laea.spatial_ref = meta['crs'].wkt
    laea.GeoTransform = " ".join(str(x) for x in meta['transform'].to_gdal())

# Second, adding x and y coordinates array to netcdf file
with netCDF4.Dataset(file2, mode='w') as src:
    # Create spatial dimensions
    # x
    x_dim = src.createDimension('x', meta['width'])
    x_var = src.createVariable('x', np.float32, ('x',))
    lon = meta2dim(meta, 'width')
    x_var[:] = lon
    # y
    y_dim = src.createDimension('y', meta['height'])
    y_var = src.createVariable('y', np.float32, ('y',))
    lat = meta2dim(meta, 'height')
    y_var[:] = lat
    # Create temporal dimension
    t_dim = src.createDimension('time', None)
    # chlor variable
    chlor = src.createVariable('chlor_a', np.float32, ('time','y','x'),
                               zlib=True)
    chlor.grid_mapping = 'laea'
    chlor[0,:,:] = array
    # laea variable to store projection information
    laea = src.createVariable('laea', 'c')
    laea.spatial_ref = meta['crs'].wkt


texts: Does anyone know the accuracy difference between the snapping used for adding a feature versus using the node tool to move a node in QGIS?

texts:
I have a point layer that I create points in first. Then I have a line layer that I create a route that snaps to each of these points using the add feature tool. When I commit to the database, I have a trigger function that creates a second line layer that by slicing the route at each of the points and stores those segments.

(The points are utility poles, the first line layer is a fiber route, and the second line layer is for creating a CAD like dimensional drawing of the distance between the poles)

In my trigger function, I use st_intersects() and it works beautifully whenever a route segment is touching at least one point that was snapped to during initial digitizing but it ignores any new points.

However, if I move a route segment around such that it is touching new poles but not ones that it was snapped to during digitizing, the segments in the second layer do not get regenerated.

Based on my testing, I suspect that the snapping accuracy between digitizing and moving a node using the node tool are different and st_intersect has an really tight tolerance.

Any idea's on how to solve my dilemma?


texts: Openlayers loop through json array and output result with as javascript

texts:
Connection.php file which stored database connection code with json:

<?php
$username = "root";
$password = "";
$database = "roadmap";
$url = "localhost";

$con = mysql_connect($url, $username , $password);
if (!$con)
{
die('Could not connect:' . mysql_error());
}

mysql_select_db($database, $con) ;
$data = array();
// query your DataBase here looking for a match to $input
$query = mysql_query("SELECT * FROM students");
while ($row = mysql_fetch_assoc($query)) {
$json = array();
$json['Id'] = $row['Id'];
$json['Lon'] = $row['Lon'];
$json['Lat'] = $row['Lat'];
$data[] = $json;
}
header("Content-type: application/json");
echo json_encode($data);
mysql_close($con);
    ?>

In separate file index.php:

$.getJSON("connection.php", function(data) {
            addMarker( data[0].Lon,data[0].Lat ,data[0].desc);
            addMarker(-0.13264,51.50918 , );
            addMarker( -0.12498,51.50807 , );
            center = bounds.getCenterLonLat();
            map.setCenter(center, map.getZoomForExtent(bounds) - 1);
            zoom = map.getZoom();

        });

this currently pulls the json data through and I have to use it for each individual assigning of a part from the database such as data[0].lon will get the longitude entry for the first record with in the database.

As I am going to be working with quite a few entry's. I want to loop through the json array and then out put all records as so:

addMarker( RECORD 1 );

addMarker( RECORD 2 );

addMarker( RECORD 3 );

Can anyone help me with this.


texts: Can ogr2ogr write data to a file geodatabase feature dataset?

texts:
I am dumping postgis tables into a file geodatabase and it works great--the result is a  file geodatabase full of featureclasses.

But is it possible to write these featureclasses into feature datasets (they are the same spatial reference):


texts: Combining ST Buffer and Intersect with Where clause

texts:
I'm trying to create a buffer around a single point in a spatial entity and finding all the points from a different entity that fall within that buffer.
I tried two approaches. The code below runs without error but returns no records, even though I know from looking at the data in ArcMap that there should records should come up (increasing the buffer to a huge number also didn't return any records). Removing the last line seems to create a buffer around every point in the incidents entity and finds every point in the hydrants entity. I want to figure out how to do this around just one point in the hydrants entity.

SELECT incidents.id_incident, hydrants.OBJECTID_1
FROM zelaya2.dbo.incidents, zelaya2.dbo.HYDRANTS
WHERE Shape.STBuffer(1000).STIntersects(incidents.geog) = 1
AND hydrants.OBJECTID_1 = 177619;

Also tried setting a temporary variable as the buffer, but again got no records returns and no error message.

DECLARE @g geography;
Set @g =
(SELECT hydrants.Shape.STBuffer(1000)
FROM zelaya2.dbo.hydrants
WHERE OBJECTID_1 = 177619)

SELECT *
FROM zelaya2.dbo.incidents
WHERE @g.STIntersects(incidents.geog) = 1

EDIT: I confirmed that both layers are SRID 4326


texts: InfoWindows in Nested Polygons (Google Maps Api 3)

texts: I have a large polygon with thick borders inside which there are several smaller polygons something like multiple city polygons inside a state border. The problem with this is that the info-windows for some of the polygons don't show up when i click the polygon.(i have set up a click listener for each smaller polygon)


texts: Software/Plugin for .mpr/.mph ( German Topo Maps, not ILWIS!)

texts:
I am looking for an alternative to GlobalMapper, regarding reading .mpr/.mph-raster ( German topo maps, not ILWIS). As far as I know, only Global Mapper can  read this raster-data properly.

Ideally I need a plugin for QGIS, which can handle/import this data-type?


texts: Extracting elevation from a TIN file with a point shapefile overlay

texts:
I have a TIN file that was built off of contour polyline shapefile (I think).  I have a point shapefile with a few points.

How can I extract elevation of the TIN file where the points occur?  I am using ArcMap with python automation.  Any free software would be welcome as well.  Perhaps a transforming of the TIN to another file type is needed.


texts: Why Do Colors Look Darker When Printed on Paper?

texts:
I made a map with a terrain-colored relief-shading in the background. I then exported them as PDF from Print Composer in QGIS. But when I print the PDF onto paper:

the shaded-relief color in the background appeared way darker than it supposed to.
the colors of roads, symbols and labels in the foreground looked all right, though..

All of the colors in the PDF itself appeared correct. Is there a way to ensure my print colors look right on paper?

Appreciate any tips/documents/links that can help me get my prints right, thanks!


texts: Change position of legend in Carto Builder

texts:
I just opened an old map of mine that I had made in CartoDB Editor in the new Carto Builder. Without doing anything to the new map, I see that the legend has been moved from the lower right-hand corner to the upper left-hand corner of the map and is now covering the map description text. I am still a noob at CSS/HTML, so I'm not sure how to move my legend back to the lower right-hand corner, and I'm not finding anything online that is working. A while back I edited the legend labels, using the following HTML code. I'm assuming the legend position has to be included in this HTML code as well? Maybe as a new div class?

 Newark (EWR) Departures

 Newark (EWR) Arrivals

 LaGuardia (LGA) Departures

 LaGuardia (LGA) Arrivals

 Kennedy (JFK) Departures

 Kennedy (JFK) Arrivals

Here's my map for reference: Flights over Queens!
(Which clearly has way more wrong with it now than just the legend. I'm assuming this is all from switching over to Builder? Oy)


texts:
When importing autocad maps, sometimes topology will be chaos, this need to be amended.
And when I add new topology in QGIS, how to save that to the postgreSQL?
Is these functions included in QGIS or secondary development will be needed?


texts: Service not available: /srv/eng/rdf.search

texts:
After upgrading to GeoNetwork 4.0.1 the rdf.search endpoint is not working. Has it been removed or has something changed?
The documentation, https://geonetwork-opensource.org/manuals/4.0.x/en/api/rdf-dcat.html, has not changed.
Error message:
Service not available
The service "/catalog/srv/eng/rdf.search" does not exist or you don't have privileges to access it. Return to the search page."
catalog is our GeoNetwork instance.


texts:
I am trying to update some code from VB6 to VB.Net 2008 and have been generally successful. However, I came across an issue where, as it processes each polygon in the shapefile, it implements the following code.

Dim lngExtRing as Integer
lngExtRing = pPolygon.ExteriorRingCount - 1
ReDim pExtRings(pPolygon.ExteriorRingCount - 1)
pPolygon.QueryExteriorRingsEx(pPolygon.ExteriorRingCount, pExtRings(0))
pExtRings(lngExtRing)

Once that code runs and you hover over the pExtRings portion of pExtRings(lngExtRing) for some polygons (not all), the first popup lists a ComObject but if you expand the plus, it indicates "Nothing".

I have run this on a shapefile that has 5 polygons and the two polygons that it fails on have negative areas. There are two ways that might resolve this, but I don't know how to implement either. One was to reverse the orientation if the area was negative, but I tried that on a different shapefile and it still failed on a polygon that did not have a negative area. The other, perhaps better, approach would be to implement the RepairGeometry class of ESRI.ArcGIS.DataManagementTools, but I have no clue on how to make it happen.

If anyone has any insight or examples of how to use RepairGeometry, I would certainly appreciate being enlightened.

Thanks

Dale


texts: Allow viewer to embed OpenLayers in their own site?

texts:
I am trying to create a way for someone viewing my site to embed the map they are viewing into their own site. Not sure where to begin on this.

Any help?


texts: Displaying X, Y coordinates on ArcGIS not working

texts:
I am trying to geocode some coordinates I have from my Excel sheet onto ArcGIS but everytime I do so, it doesn't show up on the screen. When I right click on the sheet layer to see where it went, it shows my points as this:

Here is also a photo of what my attribute table looks like:

I have tried multiple times removing it, restarting ArcMap and  but this still happens. Before, it used to display all my longitude points as 0, but since then I have tried tweaking around on my Excel sheet.

This has never happened before. Could it have something to do with ArcMap not being able to read the decimal placements since they are not aligned perfectly.

I also made sure to pick my X value to be my longitude, and Y value as latitude.

Can anyone explain why this is happening?


texts:
Is it possible to set different geographic permission to different users/roles in Geoserver?

I guess one solution would be to create different layers where each layer only shows information in that specific region. And then you set the permission on the layer.
But if you have many different "geographic permissions" that would take a lot of time to set up.

So is there a way, in GeoServer, to set a geographic boundary for a user to only view information within that boundary?


texts: Graphic style not working for line feature in GeoServer?

texts:
I am trying to use a graphic style (in ysld format) to style a line feature. The end goal is to style the line as dashed line with a picture or a circle/times in the gaps of the dashes. I tried to use the code from the ysld cookbook example but it does not work. It shows the dashes as it should but the gaps do not have the picture symbol or even a circle/times symbol. Instead it just shows that as a small dash with default black color.

Is there something I am not doing right?

Here is the sample code I am using

title: 'YSLD Cook Book: Dash/Symbol line'
feature-styles:
- name: name
  rules:
  - symbolizers:
    - line:
        stroke-color: '#0000FF'
        stroke-width: 1
        stroke-dasharray: 10 10
    - line:
        stroke-color: '#000033'
        stroke-width: 1
        stroke-dasharray: 5 15
        stroke-dashoffset: 7.5
        stroke-graphic-stroke:
          size: 5
          symbols:
          - mark:
              shape: circle
              stroke-color: '#000033'
              stroke-width: 1

The image above should be the result but instead I just get a black dash instead of the circles in the gaps.


texts: Automatically labeling layers according to raster characteristics

texts:
The question is oriented more towards end user and came up during the preparation of field maps, which have to show previously mapped biotopes and their ID on an aerial image. The map is used to rework the old polygon layers according to the current situation in the field.
The problem is that with any automatic labelling method used in the QGIS label toolbox, the labels still hide relevant areas on the printed map which are needed for orientation or for spotting potentially relevant areas or changes of vegetation patterns.
I wonder whether there is a possibility to provide rules to QGIS to place labels of a polygon layer at indicated areas according to a raster image.
The relevant background would be an aerial image, used at a 1:5000 scale in the printed map.
The labels have to be applied according to a couple of tricky rules:

the labels should not overlay distinct features in the aerial photo, e.g. single trees, small colourshifts in grassland that indicate change in vegetation pattern, small bushes, rivers, small areas with open water, ...

they are allowed to be placed on e.g. forests, housed areas, fields, so to say either very homogeneous regions or assigned areas which are of no interest

they should not overlap with each other

they should be as close to the polygon as possible

if the aerial photo indicates homogeneity and the polygon is big enough to incorporate the label, the label should be placed inside the polygon

I have some programming skills in R and could get help with python, too.
But I have no clue whether that is possible at all.
My first idea was to try to classify the rasters based on areas where labeling is possible and where it is not. This would need a bit of a training of an algorithm I guess, as the rules are probably a bit implicit.
But this can only help if labelling according to another dataset, in this case a raster, is possible at all.


texts: Why does leaflet treat an svg file differently in an image overlay than a png file and how can I prevent this?

texts:
I have a vector graphic in adobe illustrator format (ai).

I would like to partially (only some layers) display this graphic as a basic map. I exported most of the layers of this vector graphic as PNG-file and created tiles using Maptiler. I cannot display this graphic in an image overlay because it is too large. That's why I chose the detour and converted the graphics to tiles.

I saved the other layers as SVG-files and show them as an image overlay.

The PNG and the SVG-files have the same size ratio and the same bounding box. The PNG is 9726px  7950px and the SVG 4863,5px x 3975.7 px.

<svg version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px"
 width="4863.5px" height="3975.7px" viewBox="0 0 4863.5 3975.7" style="enable-background:new 0 0 4863.5 3975.7;"
 xml:space="preserve">

Unfortunately, the SVG files do not fit exactly over the tiles I created. There are small shifts. You can see it in the picture.

If I save the SVG-File as an PNG and load the PNG in the image overlay everything is fine!

I think there is something wrong with the projection. When creating the tiles, I tested with different coordinate systems. Among other things with EPSG:4326 (Bounding Box 7.088 51.401 7.365 51.541) and with EPSG:3857 (Bounding Box 789032.55 6692534.87 819868.05 6717554.10).

This is my code:

You can see a temporarily live demo here: https://fahrplan.astrid-guenther.de/

I asked this question before here: https://github.com/Leaflet/Leaflet/issues/7051


texts: ArcGIS 10.2 Combine Two Feature Classes in One File Geodatabase with Attachment (Pictures)

texts: I am working in ArcGIS 10.2, and have point data collected from 2 separate tablets with the same database schema.  Both of these files have attachments (pictures) with them.  I am trying to combine these two files and can't seem to find a way to keep the attachments when I do combine them.  I have tried creating a new file geodatabase and exporting the files as XML's.  I have tried loading the data into a new file geodatabase and nothing seems to work. Has anyone run into this, or know a way to combine these two feature classes?


texts: Where to store geodata to be layered over an OSM tile map?

texts:
I'm working on a project for which most data can come straight from OSM. But some data is not relevant to OSM and will need to be stored separately then layered over the top. Probably the tiles will be generated in TileMill, and served with Leaflet. Either I'll do the layering within TileMill, or on the fly in Leaflet.

The question is - where to store the extra geodata, and how to manage it? The data will mostly be POIs, some of which will have lengthy explanatory text and even photos associated with it. As a ballpark, let's say there will be 1000 POIs in total. (It's not necessary that the photos be stored in the same system.)

This data will need to be maintained by a loose team of volunteers, so I'm looking for some kind of (free) platform that supports online editing. It would be great if it were possible to automatically layer a whole set of related geo objects, rather than having to specify them all individually.

I'd also like to avoid having to run a server of my own just for this purpose.

One possible solution I thought of was using Google Maps's "my places" (or whatever they call it now) to create the objects, then download them say to Dropbox as .kml's. That seems a bit clunky. What would be some better ways?


texts:
I have a vector symbol that I am using rule based style with. In the Rule properties tab I click on the SVG marker and change the "rotation" which is a field I have made to handle to the rotation angle I want per each row for this feature. Am I setting this up wrong?

Every record has a slightly different Rotation so when I look at the symbols they should all not be straight up and down.


texts: Elevation profile geoprocessing service - error 000840

texts:
I am trying to create an Elevation Profile service using the Python code found on the Setting up an on-premise elevation Profile service page:

And I followed all the instructions. I am able to run the tool with no errors in ArcGIS and I can publish the geoprocessing service, but when I run the gp service inside ArcGIS I get an error saying:

Error executing tool.: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a TIN Layer.
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Terrain Layer.
ERROR 000840: The value is not a LAS Dataset Layer. Failed to execute (InterpolateShape). Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a TIN Layer.
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Terrain Layer.
ERROR 000840: The value is not a LAS Dataset Layer. Failed to execute (InterpolateShape). Failed to execute (Profile). Failed to execute (Profile).

The problem seems to be in the Interpolate Shape analysis in line 307 in the Profile Tool.py but I am not sure what is wrong with it. Do you guys have any ideas?


texts:
is there a way in ArcGIS to save ModelBuilder (with all components) so that later on I can work on it?

There is 'save' in ModelBuilder but when I wan to import it to ModelBuilder it gives me all as one component instead of all original components.


texts: PostGIS Geometry type not matching column type in QGIS

texts:
I'm trying to import some polygon shapefiles to Postgis in QGIS via SPIT. One of them can't be imported and returns this error:

ERROR: Geometry type (Polygon) does not match column type (MultiPolygon)

I've tried to change attribute "Feature Class" in that table in SPIT widnow, from MULTIPOLYGON to POLYGON, but nothing happened.
Is there any way to make PostGIS to accept both types (polygon and multipolygon) or can I convert shapefile geometry from polygon to multipolygon?
I've tried QGIS 2.0.1 in Windows and QGIS 2.3 in ArchLinux.


texts: Cost Attributes- Network analyst for Arc GIS 10.2

texts:
Can anyone tell me whether Network analyst can consider multiple cost attributes? I am looking for most ecologic route (not shortest), therefore I have assigned a value of pollution to each road segment. I have 2 layers in my analysis road and cycle layer, both have these pollution values for each segment (cycle routes have mostly lower values) and I selected this pollution attribute as a default cost attribute but would also like to consider another attribute in the best route analysis.

Can this be done or does network analyst always consider just one cost attribute?


texts: Is it possible to realign polylines from KMZ in Google Earth?

texts:
I have a street network that has been extracted from AutoCAD DWG's and converted to a KMZ file. The resulting polylines, when launched in Google Earth, do not visually align with the imagery.

I am not sure what the source data is, and not really concerned with anything other than realignment (transformation?). I have ArcExplorer but no other GIS at the moment other than AutoCAD and Microstation, which aren't true GIS platforms.

Apologies for the rambling, post, thanks.


texts: Why does writing large shapefiles crash QGIS?

texts:
I'm trying to write large shapefiles with QGIS 2.0 API, using Python, which crashes QGIS. The data is composed of simple two point lines and a data table of five columns. A thousand features or so works fine, but at 10k the program crashes for sure. Perhaps some sort of writing in chunks is a solution, or the problem lies in my code :

def write_lines (file_name, data_list, coordinate_ref_system):

    fields = QgsFields()
    fields.append(QgsField("Data_1", QVariant.String ))
    fields.append(QgsField("Data_2", QVariant.String))

    writer = QgsVectorFileWriter( file_name + ".shp", "CP1250", fields,
                          QGis.WKBLineString, coordinate_ref_system)
    if writer.hasError() != QgsVectorFileWriter.NoError:
        QMessageBox.information(None, "ERROR!", "Cannot write file (?)")
        return 0

    for r in data_list:
        # create a new feature
        feat = QgsFeature()

        l_start=QgsPoint(r[0],r[1])
        l_end =QgsPoint(r[2],r[3])
        feat.setGeometry(QgsGeometry.fromPolyline([l_start, l_end]))

        feat.setFields(fields)
        #feat.setAttributes([ str(r[4]), str(r[5]) ]) DOESN'T WORK?

        feat['Data_1'] = r[4]
        feat['Data_2'] = r[5]

        writer.addFeature(feat)
    del writer

Thanks!


texts: Generating a routing database for the entire world while keeping the tags

texts:
I've been looking into this matter for a few days now, but I haven't come up with a working solution.

Basically, what we are trying to do is create an entire world database, preferably in PostgreSQL, using PostGIS. Our world file is the whole planet data, in a pbf format (~34GB).

So far, I have tried using the osm2pgrouting and osm2po tools to generate this database, but both tools failed for the setup that we are using.

Our machine is a linux x64 ubuntu-server, with 20 GB of RAM plus 20 GB of swap area. We are using an HHD of 2TB for this.

In my last try with osm2po, I tried running it with 18 GB of ram (-Xmx18432m), but it failed after some time. I will be trying again with more memory allocated to the java virtual machine and report the results.

But mainly, my main questions here are:

Is it possible to generate a routing database for the entire world with these tools and our current machine?
Is RAM allocation the only optimization that I can do to allow osm2po to handle large files?
Can osm2po save nodes tags in the PostgreSQL database? I know this is possible for the ways table, since it already exists in the database and even has a column called osm_meta, and I also know that it is possible to retrieve this info while reading osm data file, but my guess is that this data will only be available to the front-end part of osm2po.


texts: Does each road segment need to be digitized into a separate entry for a road network?

texts:
I am trying to understand this code in pgRouting for creating topology:

CREATE OR REPLACE VIEW road_ext AS
   SELECT *, startpoint(the_geom), endpoint(the_geom)
   FROM road;
CREATE TABLE node AS
   SELECT row_number() OVER (ORDER BY foo.p)::integer AS id,
      foo.p AS the_geom
   FROM (
      SELECT DISTINCT road_ext.startpoint AS p FROM road_ext
      UNION
      SELECT DISTINCT road_ext.endpoint AS p FROM road_ext
        ) foo
   GROUP BY foo.p;

CREATE TABLE network AS
   SELECT a.*, b.id as start_id, c.id as end_id
   FROM road_ext AS a
   JOIN node AS b ON a.startpoint = b.the_geom
   JOIN node AS c ON a.endpoint = c.the_geom;

From this code, can anyone explain this code a bit please. I am confused and want to know if it generates entries of individual road segments in separate rows even if I digitize some road segments into a single road.

Clearly, I want to know what way we need to digitise roads for creating a road network? Shall we need to split a long road into new segments at  every connection to a small road or do the below code generates it automatically? Please help!!

As the code only says startpoint and endpoint, I dont get how it can identify the intermediate nodes where there are connecting roads.

For Example, in the below image

I digitised the long road in a single entry. Now if I run the above code, will it generates the individual line segments?


texts: Can i combine DXF, Google Earth Models, Shape Files and Rasters in QGIS?

texts:
I'm an architect, trying to bring a number of data sets together for use in a CAD package (AutoCAD).  I have some georeferenced images (.jpg, .jgw.), Some model files downloaded from GoogleEarth (.kml and .skp), as well as digimap data such as vector maps (.shp) and contour profiles (.dxf). AutoCAD does not support georeferencing, so i'd like to combine this information in another package, such as QGIS, and then export it for working on within AutoCAD.

I've been playing with this, and not had much look.  I can't see how to import model files (.skp) or lines (.dxf), and i'm not sure if this is possible.  Also, export is failing for the DXF export of vector lines from things like the OS local map.

Any ideas about these problems, or other suggestions as to programmes that could bring these various kinds of information together?


texts: Querying contour lines for elevation at point using PostGIS?

texts:
Is there a way to interpolate between these contour lines for a given point using postgis?

I've come up with the following query to gather the potential contour lines of interest:

        $sql = "SELECT
                    osm_id,
                    ele,
                    ST_AsGeoJSON(way) AS geojson,
                    ST_distance(way, ST_Transform(GeomFromText('POINT(".$lon." ".$lat.")',4326),900913)) AS dist,
                    ST_ClosestPoint(way,ST_Transform(GeomFromText('POINT(".$lon." ".$lat.")',4326),900913)) as closest,
                    ST_Azimuth(ST_Transform(GeomFromText('POINT(".$lon." ".$lat.")',4326),900913), ST_ClosestPoint(way,ST_Transform(GeomFromText('POINT(".$lon." ".$lat.")',4326),900913)))/(2*pi())*360 as angle
                FROM
                    planet_osm_line
                WHERE
                    way && ST_Transform(GeomFromText('POINT(".$lon." ".$lat.")',4326),900913)
                ORDER BY
                    dist
                LIMIT 10";

The post processing of these results I will just describe rather than posting my convoluted code.  It seems the real crux of the problem is finding the two contours of interest.  Using the two closest contours isn't always the best or even accurate solution as they may both be on one side of the point of interest.

What I did was accept the closest contour as a keeper. Then check the azimuth between the point of interest and the closest point and compare this azimuth to the similiar azimuth calculated from the first accepted contour.  The closer the delta between these two was to 180 the higher the probability that this was my second contour.


texts: Sentinel-1 radar data all zeros with Python SentinelHub API

texts:
I am requesting Sentinel-1 SAR data but rather than throwing an error the code is just returning zeros.
My request:
bbox = BBox([-145.8, 60.96, -146.5, 61.38], crs=CRS.WGS84)
bbox_size = bbox_to_dimensions(bbox, resolution=100)
time_interval = ('2018-07-01', '2018-07-16')

evalscript_all_bands = """
    //VERSION=3
    function setup() {
        return {
            input: [{
                bands: ["HV"],
            }],
            output: {
                bands: 1,
            }
        };
    }

    function evaluatePixel(sample) {
        return [sample.HV];
    }
"""

request_all_bands = SentinelHubRequest(
    evalscript=evalscript_all_bands,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL1_EW,
            time_interval=time_interval,
    )],
    responses=[
        SentinelHubRequest.output_response('default', MimeType.TIFF)
    ],
    bbox=bbox,
    size=bbox_size,
    config=config
)

data = request_all_bands.get_data()

print(len(data))

print(data[0][:,:].shape)

plt.imshow(data[-1][::],aspect='auto')

The output of this is an array of zeros - no error message. Can anybody suggest a time period or region where I can find some data?


texts: Create line segments at point coordinates in QGIS

texts:
I have a point vector layer with angle field. How can I create line segments at these points with certain length and angle?

I did this with style options but I need this physically as lines vector layer:


texts: GEMI vegetation index calculation in ENVI

texts:
I want to calculate gemi index, but I can't write its formula in ENVI, this doesn't work.

gemi formula is:

(n(1−0.25n)−RED−0.1251−RED)
n = ( 2 * ( NIR ^2 - RED ^2) + 1.5 * NIR + 0.5 * RED ) / ( NIR + RED + 0.5 )

I bring these formulas in one, but it doesn't work:

(2*(b4^2-b3^2)+1.5*b4+0.5*b3)/(b4+b3+0.5)*(1-0.25*(2*(b4^2-b3^2)+1.5*b4+0.5*b3)/(b4+b3+0.5))-((b3-0.125)/(1-b3)

In what way I can calculate GEMI?


texts:
I know this question (or a version of it) has probably been asked but I cannot seem to figure out a solution to my problem.

The core spatial operation is finding which of my points are inside polygons. I am using my own coordinates and Census Bureau shapefiles in my database.

I have a PostGIS database set up with a schema used for doing my stuff. It looks like this:

create schema temp2;

/* load data from csv into this table
*/
create table temp2.coords (
id integer primary key,
lat float,
lng float
);

/* contains the geoms
*/
create table temp2.geoms (
id integer primary key,
geom geometry(point, 4269)
);

/* my spatial index
*/
create index temp2_geoms_idx on temp2.geoms using gist(geom);

/* save the result
*/
create table temp2.result (
id integer primary key,
countycode varchar(64)
);

The main issue is that when I run the query, it takes forever for a small amount of points. It took about 30 min for 100k points. From what I have read, this should happen way faster, and I am sure it has to do with the query that I have.

I have experimented with different spatial indices, different kinds of queries, and read the manual. I think I am misunderstanding something about how the indices work.

This is my query:

select
main.id,
c.cntyidfp
from temp2.geoms as main
left outer join tiger.county as c
on st_intersects(c.the_geom, main.geom);

Explain from the query:

"Gather  (cost=1000.00..1870856.75 rows=6683599 width=10)"
"  Workers Planned: 2"
"  ->  Nested Loop Left Join  (cost=0.00..1201496.85 rows=2784833 width=10)"
"        ->  Parallel Seq Scan on geoms main  (cost=0.00..77500.33 rows=2583333 width=36)"
"        ->  Append  (cost=0.00..0.42 rows=2 width=39218)"
"              ->  Seq Scan on county c  (cost=0.00..0.00 rows=1 width=41)"
"              ->  Index Scan using tiger_data_county_the_geom_gist on county_all c_1  (cost=0.15..0.42 rows=1 width=39230)"
"                    Index Cond: (the_geom && main.geom)"
"                    Filter: _st_intersects(the_geom, main.geom)"
"                    Filter: ((the_geom && main.geom) AND _st_intersects(the_geom, main.geom))"

Shouldn't a simple query like the one I have written be fast?

I thought the whole point of the spatial index was to speed up the "st_" functions.

Shouldn't setting the spatial index to the geometry speed this up?

How can the query be written to avoid this?

EDIT: Query was changed to:

select
m.id,
c.statefp
from tiger.state as c
inner join temp2.geoms2 as m
on st_intersects(m.geom, c.the_geom);

And was able to get a significant speed increase. About 6.2 million points in under 1 minute.


texts: Using a fixed palette range to plot an NDVI raster in R

texts:
I would like to plot an NDVI raster with rasterVis packages using the normalized colors from red to green. I don't understand how fix the color scale from 0 to 1 even if my raster value are between 0.15 and 0.40

I have tried :

library(rasterVis)

n <- 1000
col.seq <- seq(0.0001, 1, length.out=n)
brks <- c(0, col.seq, Inf)
cuts <- cut(as.matrix(stackmean), breaks = brks)
colors <- colorRampPalette(c("red", "yellow", "lightgreen"))(length(levels(cuts))-1)

levelplot(stackmean, col.regions=colors)

My scale clearly doesn't work !


texts: How to adjust/convert/transform set of points with unknown crs using reference points

texts:
I have a set of points in an unknown crs.
How can I convert (not guess) this to a given crs?
I know, or rather I can determine the actual coordinates of each point in order to have some reference points to match against.
How can I convert the rest of the points using these reference points?
I have done this with raster images (using qgis + georeferencer plugin) , but I don't seem to find any tools for doing this with points.
I don't need great accuracy.

I have managed to center it and scale it in a very naive and linear based on a single reference point (near the center) but the accuracy is pretty bad, specially with points far away from the reference point.


texts:
I did setup sucessfully a OSM Website following tutorial at https://github.com/openstreetmap/openstreetmap-website

I have now an OSM server empty of vectors.
Now I want to import some data of shp format distributed along the whole country of Brazil.
I tried using JOSM, but when I try to retrieve the data from my server it says that the area is too big (even a medium city it says that its too big.

Is there a way of importing and deleting data of large areas using JOSM?


texts: Catching shp2psql errors from python's subprocess module

texts:
I'm importing data into postgresql with shp2psql and pyhton's subprocess module.

command = 'shp2pgsql -s {SRID} -a -g geom -W "UTF-8" {SHP_path} {table_name} -I | ' \
          'psql -h {hostname} -p {port} -d {db_name} -U {username}'\
    .format(SRID=4326, SHP_path=shp_file, table_name=table_name,
            hostname=hostname, port=port, db_name=db_name, username=username)

subprocess.call(command, shell=True)

The problem is that I cannot intercept import errors like

ERROR:
current transaction is aborted, commands ignored until end of transaction block

They are happening internally and don't stop the subprocess from execution so I have no means to control the result of execution. Is there any way to raise such errors in the main code?


texts: Auto populating XY fields when new point feature digitized using ArcGIS Desktop?

texts:
Is there a way for ArcGIS Desktop to populate a field automatically when a new point is dropped on the map?

We have our map in State Plane Feet, and we want to the database to automatically add WGS XY coordinates to fields in the table when a new point is added.

Something like a calculated field, but always there?


texts: Joining Excel file with CityGML file using FME?

texts:
I want to join a excel file with a CityGML file with FME.

i want to join the excel file with the geometry of the CityGML file.

all attributes should join this file

which settings should i choose in the feature merger to get it?

in the end i want to put it all together in a ArcGIS geodatabase


texts: How to change raster cell size WITHOUT resampling? (ArcGIS Pro)

texts:
I have hundreds of rasters, all with different rectangular cell sizes. I need to convert the pixels to be square, since I will be using a variety of Spatial Analyst tools such as Raster Calculator. However, I do not want to use the Resample tool as I do not want to reduce data accuracy. I am okay with the extent/size of the raster changing. How would I do this in ArcGIS Pro?

Edited for more info:
although my rasters are all different sizes, when converting them to square cell sizes, it is ok for them to continue to be different sizes (just as long as they are square). Basically, I have hundreds of rasters from different areas, and my workflow is:

Multiply classRasterA1 by certaintyRasterA1 (Area 1)

Multiply classRasterA2 by certaintyRasterA2 (Area 2)

and so on.


texts: Computing producer and user accuracy for each class in Google Earth Engine supervised classification

texts:
I know I can get a confusion matrix and an overall accuracy, but is there a way to calculate producer and user accuracy for each class separately from the Google Earth Engine?

This is my validation dataset: https://code.earthengine.google.com/?asset=users/RohitNandakumar/accuracy

This is my training dataset: https://code.earthengine.google.com/?asset=users/RohitNandakumar/train

var sent = ee.ImageCollection("COPERNICUS/S2_SR"),
    table = ee.FeatureCollection("users/RohitNandakumar/train"),
    geometry =
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[79.34242573948666, 29.562958506122595],
          [79.34242573948666, 29.296831501244323],
          [79.79217854710384, 29.296831501244323],
          [79.79217854710384, 29.562958506122595]]], null, false),
    accuracy = ee.FeatureCollection("users/RohitNandakumar/accuracy");
var land = sent.filterDate('2019-04-1','2019-11-01')
.filterMetadata("MGRS_TILE","equals",'44RLT')
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE",5))
.sort('system:time_start');
var image = ee.Image(land.toList(land.size()).get(1)).clip(geometry)//.aside(print);

//supervised classification
var label = 'MC_ID';
var bands = ['B2','B3','B4','B5','B6','B7','B8','B8A','B9','B11','B12']
var training = image.select(bands).sampleRegions({
  collection: table,
  properties: [label],
  scale: 10
});
var training = training.randomColumn().filter(ee.Filter.lt('random', 0.1))
var trained = ee.Classifier.cart().train(training, label, bands);
var classified = image.select(bands).classify(trained);

var validation = image.select(bands).sampleRegions({
  collection: accuracy,
  properties: [label],
});
var validated = validation.classify(trained);
var testAccuracy = validated.errorMatrix('MC_ID', 'classification');
print('Validation error matrix: ', testAccuracy);
print('Validation overall accuracy: ', testAccuracy.accuracy());


texts: Stack 3 layers and display them as one in GEE

texts:
I have used this code https://www.gears-lab.com/intro_rs_lab8/ in order to learn and practice my SAR-1 skills.
I have created an almost similar code, but few things are not clear to me regard this code.
In the end of the code, I am supposed to display in RGB 3 different layers. How does it work? I get only one image, so what is the meaning of display each one of them as RGB in one image?

When I use the inspector in order to check the values I gave in the RGB stack (seasonal stack) I have "layer" that not suppose to be there- VV. That confuses me even more- why don't I get the 3 layers I stack together? (vv1,vv2, and vv3)

I would say this relates to the first question. of how it displays 3 images as 3 different band in one image.

This is the code I used, it has error with the VH for a reason I don't know yet

var geometry=geometry2;
// Filter the collection for the VV product from the descending track
var collectionVV = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterBounds(geometry)
    .select(['VV']);
print(collectionVV);

// Filter the collection for the VH product from the descending track
var collectionVH = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterBounds(geometry)
    .select(['VH']);
print(collectionVH);

Map.centerObject(geometry);

var clippedVV=collectionVV.map(function(im){
   return im.clip(geometry);
});

var clippedVH=collectionVH.map(function(im){
   return im.clip(geometry);
});

var VV = clippedVV.median();
var VH= clippedVH.median();

// Adding the VV layer to the map
//Map.addLayer(VV, {min: -30, max: -1}, 'VV');
//Map.addLayer(VH,{min:-30, max:-1}, 'VH');

// Create a 3 band stack by selecting from different periods (months) for VV

var VV1 = ee.Image(clippedVV.filterDate('2017-10-01', '2017-11-30').median());
var VV2 = ee.Image(clippedVV.filterDate('2017-12-01', '2018-12-31').median());
var VV3 = ee.Image(clippedVV.filterDate('2019-01-01', '2019-02-28').median());

// Create a 3 band stack by selecting from different periods (months) for VH

var VH1 = ee.Image(clippedVH.filterDate('2017-10-01', '2017-11-30').median());
var VH2 = ee.Image(clippedVH.filterDate('2017-12-01', '2018-12-31').median());
var VH3 = ee.Image(clippedVH.filterDate('2019-01-01', '2018-02-28').median());
//Add to map
Map.addLayer(VV1.addBands(VV2).addBands(VV3), {min: -30, max: -1}, 'Season compositeVV');
Map.addLayer(VH1.addBands(VH2).addBands(VH3), {min: -30, max: -1}, 'Season compositeVH');

Map.addLayer(VV1,{min:-30,max:-1},'Season composite inly vv1');
Map.addLayer(VV2,{min:-30,max:-1},'Season composite inly vv2');
Map.addLayer(VV3,{min:-30,max:-1},'Season composite inly vv3');


texts: Highlighting or marking property using node id in leaflet?

texts:
I am working on openstreetmap using leaflet.js, I found some difficulties in highlighting a way(road) using it's node ID of OSM data.

I need to know is it possible or not?


texts: Offline Python or C++ library for plotting locations on map?

texts:
I'd like to plot locations using longitude and latitude coordinates using Python or C++, is there a good library for such a task?

I know that matplotlib is an option, I'm wondering if there are any others?


texts:
I was just asked in a interview what a dangling nodes was.  I had no idea and asked a former co-worker.

I was informed that a significant amount of time was spent fixing these.  I was confused how this could become a problem and came up with few quick solutions.

Go through the system and identify all endpoints.
Reduce the list of endpoints to endpoints that are terminals. (optional just less work)
Eliminate all endpoints that aren't within a range of another object (within the same feature type)
Create a report
Have a system to mark the entry as resolved, if it is a valid entry (so it is not reproduced in the report).  For example two roads that are really close but a 1 meter ditch separates them)
Repeat regularly.

Solution 2:

Check for system percolation
Mark all entries that do not percolate within a small area
Steps 4 through 6

This would produce a larger amount of false positives, but it's running time is in order (N log N)

We were just going through the map block by block looking for them.  I must be missing something since I can think this problem can be solved relatively simply.  I realize that the initial round would be quite a task, but with a system in place maintenance and monitoring would be minimal.

What am I missing that makes dangling nodes a large issue?


texts: Which polygons surround this point? A practical case of political geography

texts:
I'm starting a project from scratch but am budget constrained. In plain language I've got no money. What I want is to provide an address and return what governments have jurisdiction there. The identity of the governments will be coming from the US Census Bureau, though I am not sure they have borders of these governments.

I'm currently considering converting the address to lat/lon coordinates and defining that as a point and defining all ~90k US governments as polygons and then running a search. The large majority of governments are going to be in a different state so the ~90k number will be cut down by a lot just by filtering by state. Is setting up 50 layers (or 50 database tables) with all governments present in a state and then running the analysis only in that layer realistic? Or am I going to need to make enough layers to do it by counties? I'm trying to get at how much calculation am I avoiding if I restrict the geography of the polygons. How much restriction would likely be optimal?

When the GIS database is set up, people will be able to query for "their" results, more exactly the results of the address they put in. I'm also going to create a background task to process addresses ahead of time. The more efficient the calculations, the more I can precalculate.


texts:
I'm using QGIS 2.18.14 and I have currently got a map tip in a shapefile of freehold land titles that looks like this:

Both, the title plan and register for the land title have been referenced to a hyperlink on dropbox.

The attribute table for the shapefile is below:

My CSS / HTML for the map tip looks like the following:

 [% "Title No" %]

Proprietor

 [%"Proprietor"%]

 Area

 [%"Ac" || ' Ac'%]

 Title Plan

Click Me

 Register

Click Me

You will notice I have not got a link to the title plan ending '11'. I will be uploading this QGIS project to QGIS Cloud and want to make it as professional looking as possible. How would I go about omitting the row "Title Plan" in the map tip represented below when the entry = 'NULL'?

 Title Plan

Click Me


texts:
In creating an export from QGIS 3.8.0 Zanzibar to Qfield using the sync plugin, I am encountering a strange error.

Each time I run the export, several of the layers I have in a geopackage in QGIS file do not correctly copy over to the newly created *_qfield.qgs file.

The odd thing is that I can't systematically reproduce the error; each time I test the export, I have different layers that successfully export or fail. There are ten layers in the geopackage.

Does anyone have thoughts on how to troubleshoot this?

I am pretty new to QField.


texts:
I would like to streamline the render of a polygon grid vector layer as a simple line (no fill) via Python coding. I am I bit confused by the implementation referenced by the API.

Is it possible to resolve some simple examples?

Here is my first blind go at implementation (currently non-functioning) -

symbol_gridLayer = layer.rendererV2().symbols()[0].changeSymbolLayer(0, QgsSimpleLineSymbolLayerV2(black, 2)


texts: Do Postgresql sequences work with SDE versioning?

texts: I am migrating data from an older GIS system into an SDE system with PostgreSQL as the database engine.  The current system generates a unique tag for several layers, I want to re-create this behavior using sequences in the PostgreSQL database.  I have it working in a single user environment, but will it work in a versioned environment?  Right now the field is populated by using Postgres' nextval() function as the field's default value.  When the tables are versioned do those defaults go with them?


texts: The first layer does not appear on the map using WMS in openlayers. What I'm missing?

texts:
I'm trying to shown on the map the layer with the url  , but I can not find the solution. However with the second layer it works. What I miss? And why?

I made this fiddle:

And yes, I have tried to change the projection, and even remove it.The layer would have to appear to the northwest of Africa, Canary Islands

EDIT:
My bad sorry

link: http://idecan1.grafcan.es/ServicioWMS/MTL

fiddle: https://jsfiddle.net/nzx9fvz5/2/

var projec = new ol.proj.Projection({
  code: "EPSG:32628"
});
  var map = new ol.Map({
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      }),
      new ol.layer.Tile({
        source: new ol.source.TileWMS({
          url: "http://idecan1.grafcan.es/ServicioWMS/MTL",
          params: {'LAYERS': "LIDAR_MTL"},
          projection: projec
        })
      }),
      new ol.layer.Tile({
        source: new ol.source.TileWMS({
          url: "https://mesonet.agron.iastate.edu/cgi-bin/wms/us/mrms.cgi",
          params: {'LAYERS': "mrms"}
        })
      })
    ],
    target: 'map',
    view: new ol.View({
      center: [0, 0],
      zoom: 2
    })
  });


texts: Optional Parameter becoming required in ArcPy?

texts:
I would like to make a python script (.pyt) where i have 2 input parameters:

If the first parameter (which is a boolean) is set TRUE, therefore it is checked, the 2nd parameter should become required (if not, the 2nd parameter should be optional).

This is the used code so far:

def getParameterInfo(self):
    param0 = arcpy.Parameter(
    displayName="BooleanTest",
    name="in_features",
    datatype="Boolean",
    parameterType="Required",
    direction="Input")

    param1 = arcpy.Parameter(
    displayName="Test",
    name="testfield",
    datatype="Field",
    parameterType="Optional",
    direction="Input")

    params = [param0,param1]

    return params

def updateParameters(self):
    if self.params[0].value == "True":
        self.params[1].parameterType = "Required"
    else:
        self.params[1].parameterType = "Optional"

    return

what am I doing wrong? I'm new to Python.


texts: I am trying to create a raster data set from county level polygon data for the US.  I have successfully created one raster from a dataset, but when attempting another, it does not work.  I am wondering if I am missing something.  It's in ArcMap 10, I have tried moving the cell size up to 100 (it took the recommended cell size for the other data set) and have also changed any NULL values to zero. I get the "999999 : Error executing function".  I have also tried the Feature to Raster tool.  Same error message.  Any ideas or troubleshooting tips would be great!  Thanks.


texts: How to interpret Tasseled Cap Greenness for Landsat 5?

texts: I have amassed a few promising literature articles, and looked at the help files for ENVI, ArcMap and PCI, and no one really goes into great depth on how to interpret spectral values for Tasseled Cap greenness. More specifically, what is the range of healthy vegetation in the greenness band? I know for NDVI it is 0.3-0.5. I'm getting different ranges of values according to different scenes for TC-green, so is the range scene-dependent?


texts: How to create an interactive map with jQuery and HTML?

texts:
I want to develop a map in JQuery and HTML. It needs to show world country borders. It must be enabled mouse hover effect and also popup (when mouse over) in which I can add my own text.

What JS plugin do you suggest? I tried JVectorMap but it has only hover effect and I couldn't make popup window with own text.


texts: Using ArcGIS Desktop to zoom to next selected feature?

texts:
How could I adapt the following code to enable me to zoom to each individual selected Geometry one after the other?

import arcpy
mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
df.zoomToSelectedFeatures()
arcpy.RefreshActiveView()


texts: Script to save lots of shapefiles into different CRS in QGIS

texts: I have a lot of shapefiles that I need to save in a different CRS... is it possible to automate this process with a script in QGIS?


texts: Lock the zoom level of WMS feeds when printing a layout in QGIS 3

texts:
This is similar to the following question which was asked 1 years ago and one asked 5 years ago.

Is there now a way to lock the zoom level of a WMS feed in a layout? It's not an XYZ layer so I can't use the strategy mentioned in the above posts.

I want it to show as below

But it always prints at a different scale level.

I have the following settings in the print to pdf setting


texts: Cannot Publish from ArcMap to ArcGIS Online?

texts:
I am attempting to publish an MXD as a map service to ArcGIS Online from ArcMap Desktop 10.3.

I have a home licence and have a Organization account on ArcGIS Online with 60 credits left.  I have published services in the past few months.

When I attempt to Publish a Service from ArcMap, I am requested to create a Connection.  (I already have one working Connection as "User" but nothing as "Publisher".)

I attempt to create a Publish connection to:
http://services.arcgisonline.com/arcgis/services

which I believe is the correct URL, however I get the error message:
"Error: Proxy Server got bad address from remote server(verify the server is running)."

Fiddler shows: [Fiddler] The connection to 'services.arcgisonline.com' failed. Error: TimedOut (0x274c). System.Net.Sockets.SocketException A connection attempt failed because the connected party did not properly respond after a period of time, or established connection failed because connected host has failed to respond 125.56.183.230:8080

I am unsure of where to go from here.  I am not using a proxy server. I have searched for this problem and found suggested solutions. I have applied the Internet Explorer fix and the ArcGIS Administrator fix but these have not worked.


texts: What C/C++ libraries exist to transform vectors between map projections

texts:
Has anybody written a C/C++ library (preferably one that is aware of Proj.4 definition strings) that has a set of routines to convert vector quantities (e.g., wind speed and direction) between map coordinate systems? That is, routines to do calculations described in something like:

http://www.dtcenter.org/met/users/docs/write_ups/velocity.pdf

I have found the EMAPF routines, but was wondering whether there was anything Proj.4 aware.

Gordon


texts:
How can i display an attribute table so it does not disappear from the screen, so i work on an open map and see the attribute table simultaneously?

When i work on a map and has also opened attribute table, the attribute table disappear when i return to work on the map.

I can minimized the map screen and work with the attribute table simultaneously,

but i prefer to work on a map with on a full screen. Is there a way to to do it?


texts: Extracting specific geometry point from MultiPoint geometry using PostGIS?

texts:
I am working with multipoint geometries in postgis and would like to extract the geometry point associated with the greatest y value of the multipoint geometry. How may I accomplish such a task?

For example, if my multipoint geometry is defined as follows:

'MULTIPOINT(2 5, 3 1, 4 0)'

Then how may I gather the geometry point associated with the greatest y-value?

The following query will tell me that the greatest y-value is 5:

SELECT ST_ASTEXT(MP), ST_YMAX(MP) FROM
(SELECT ST_GeomFromText('MULTIPOINT(2 5, 3 0.9, 4 1.2)') as MP) as foo;

However, I would like to create a query that returns 'POINT(2 5)' rather than the YMAX value of 5.

Any thoughts?


texts: Openlayers: How to make layer come to top after choosing it in layer switcher?

texts:
I have a base layer and 3 image layers . I also have setup a layerswitcher so that user can choose which layer is to be overlayed on the base map and I have 3 sliders such that user can choose the opacity of each layer. The problem as of now is that because of layer indexing, for the user to see the image layer right above the base layer other layers have to be switched off. I want it such that when any of the layers is switched on it automatically comes to the top.

This is a link to the site.

And my code is below:

        var map = new OpenLayers.Map('map');

        var NaturalEarth = new OpenLayers.Layer.WMS(
        "Satellite","http://arcims.csl.gov.uk/wmsconnector/com.esri.wms.Esrimap/naturalEarth?",
    {layers: 'naturalEarth'}
    );
        map.addLayer(NaturalEarth);
        map.zoomToMaxExtent();

        var imgLayer1 = new OpenLayers.Layer.Image(
            'T2',
            "http://wrf1.geology.um.maine.edu/bipush/OpenLayers-2.13.1/OpenLayers-2.13.1/kml/T2.png",
            new OpenLayers.Bounds(-180, -90, 180, 90),
            new OpenLayers.Size(512, 256)
            //{numZoomLevels: 3}
     );

    var imgLayer2 = new OpenLayers.Layer.Image(
            'PWTR',
            "http://wrf1.geology.um.maine.edu/bipush/OpenLayers-2.13.1/OpenLayers-2.13.1/kml/PWTR.png",
            new OpenLayers.Bounds(-180, -90, 180, 90),
            new OpenLayers.Size(512, 256)

            //{numZoomLevels: 3}
     );
     var imgLayer3 = new OpenLayers.Layer.Image(
            'MSLP',
            "http://wrf1.geology.um.maine.edu/bipush/OpenLayers-2.13.1/OpenLayers-2.13.1/kml/MSLP.png",
            new OpenLayers.Bounds(-180, -90, 180, 90),
            new OpenLayers.Size(512, 256)

            //{numZoomLevels: 3}
     );

    imgLayer1.setIsBaseLayer(false);
    imgLayer2.setIsBaseLayer(false);
    imgLayer3.setIsBaseLayer(false);
    map.addLayers([imgLayer1,imgLayer2,imgLayer3]);
    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.zoomToMaxExtent();

  </script>

        0% <input id = "slide" type = "range"
        min ="0" max="100" step = "1" value ="100"
         onchange = "Layer1OpacityUpdate(this.value)" />100% &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

        0% <input id = "slide" type = "range"
        min ="0" max="100" step = "1" value ="100"
         onchange = "Layer2OpacityUpdate(this.value)" />100% &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

        0% <input id = "slide" type = "range"
        min ="0" max="100" step = "1" value ="100"
         onchange = "Layer3OpacityUpdate(this.value)" />100%

Thanks in advance


texts:
I have a tab file with the contour lines for an area, and also a tab file with the land polygon.

Now, I'm following a tutorial written by a colleague a couple years ago, where a thematic map is created based on those same tab files to eventually make a 3D map. It fails at the step of creating the thematic map. The generated layer looks blurred and does not stop at the land polygon.

Neither me nor my colleague can replicate the results from the tutorial. We are using MapInfo 15, and a different (newer) computer. I have also tried with an older version of MapInfo (12.5). Could this be related to our graphics hardware?

Here are the steps I followed (Mapinfo 12.5):


texts: How to access Landsat bucket from Amazon?

texts:
I tried to access amazon Landsat S3 bucket with the following python code

import boto
import boto.s3.connection
from boto.s3.connection import S3Connection

access_key = <userid>
secret_key = <pwd>
conn = S3Connection(access_key, secret_key)

nasaNexBucket = conn.get_bucket("nasanex")

for key in nasaNexBucket.list():
    print key.name

and getting this down error. Can anyone help me how to access landsat s3?

  " nasaNexBucket = conn.get_bucket("nasanex")   File
  "C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\boto\s3\connection.py",
  line 502, in get_bucket
      return self.head_bucket(bucket_name, headers=headers)   File "C:\PROGRA~2\QGISWI~1\apps\Python27\lib\site-packages\boto\s3\connection.py",
  line 535, in head_bucket
      raise err S3ResponseError: S3ResponseError: 403 Forbidden "


texts: Can typical Meets and Bounds (angles and distances) drawings be accurate in QGIS?

texts:
I’ve been working with parcel boundary shape files (.shp), standard satellite imagery, and Meets and Bounds (angles and distances) site maps / legal descriptions for a few weeks.  Wondering if those can be accurately rendered in QGIS, down to about 6 inches? (yes the readily available satellite images are very fuzzy at a close range, but a high contrast wall, driveway, roof etc. can be easily seen to violate vector property boundaries.)

What level of error would be typical of that much zoom  (down to perhaps 50 x 50 feet)

More details below:

There is the Azimuth and Distance Plugin (very popular - 46k downloads!) that allows the creation of a feature from angles and distances, however I’ve found it (for QGIS v 2+) early alpha level (no documentation, UI has features not implemented, etc. and buggy to the point of not being usable in versions 2+ of QGIS.

In Azimuth and Distance Plugin in QGIS 1.8 has less features but seems to work fine. Best of all it imports and exports from a standard text file so writing to those can avoid the tedious data entry step though the UI.

You can only construct a line feature in the plugin and that is limiting. Ideally, multiple polygons per feature are needed.

Many U.S. Counties now have .shp files available of all their property boundaries. How did these .shp files actually get created in the first place, and about how accurate are they? I’d assume they digitized the historical neighborhood blueprints (Plats) on file with the County and then hand entered textual Meets and Bounds on the remaining ones, rural ones, not on a Plat ones? That must have been a ton of hand work!

Would they have gone to the scanned original deeds and created CAD drawings (DXF format?) from the text descriptions of those? Was there any automation? Often times the local survey markers stuck in the earth are hard for people who are not surveyors to pinpoint and tag with highly accurate geo coordinates.

In many counties a grid system is used of 1 mile squares, further divided into .25 mile squares. The Point of Origin for a Plat is typically one the corners on the .25 mile squares. Would a County provide pinpoint accurate geo coordinates for those .25 mile squares? Might be great to have those as a pinpoint accurate .shp file.

I’m guessing over the typical acre property (about 200 feet square) a land surveyor working from County established fixed earth markers during the 1940s-50s-60s-70s and typically 200 foot vector segments could sight in boundary pins with a radius of around 6 inches from each node? Errors in angle measurement and distances from there might accumulate, and the entire boundary could be off by a foot or more from start to finish? In the 70s I noted a large retail store in my home town managed to be off on its structural steel and footers a full 1.5 feet! It delayed the completion of construction by several months because nothing fit. Were errors in site survey work of that magnitude fairly commonplace?

Taking the old Meets and Bounds from paper to a vector drawing file has a chance of error but as long as the original survey was high quality (for it's day) I can’t see that step as introducing much error.

However, connecting the Point of Origin and any physical land marks properly to the other map layers in the GIS strikes me as potentially difficult to get correct. Not only is the “connecting” a challenge but how about projection and zooming? Particularly, all way in (to the site plan level), it strikes me as potentially quite risky? I need to admit I know little about the different CRS and how well they scale.

Traditional CAD drawing software packages and GIS systems evolved separately over decades, each to serve their unique and traditional target audiences. Now the two disciplines seem to be becoming more integrated, particularly with on the ground surveyors using commercial grade GPS.

So, the question is not really just how accurate can QGIS be at the site plan level, it's questioning if its a good idea in general to push together the formerly distinct mapping functions?


texts: How to set raster layer name inside raster::calc function

texts:
I'm using raster::calc to make some calculations and write the output to GeoTiff in just one go. But I'd like to set the layer name of the output before writing to disk. I'm trying to do this by using setNames inside the function but it's now working. How can I change layer name inside a raster::calc call that saves the file to disk?

Here is a simple reprex

library(raster)
#> Loading required package: sp

r <- raster(ncols=36, nrows=18)
r[] <- 1:ncell(r)
names(r) <- 'my_name_is'

r2 <- calc(r, function(x) {
    y <- x * 10
    y <- setNames(y, 'my_name_is')
})

names(r2)
#> [1] "layer"

Created on 2020-05-07 by the reprex package (v0.3.0)


texts:
I am importing a .txt file into ArcMap (10.4.1). It all works fine, apart from the fact that when I open the table in ArcGIS it says the fields are string, when they should be numeric [please see image below].

The original .txt file looks like this:

How can I make sure that the values are numeric, and not strings?


texts: Access-Control-Allow-Origin in geoserver and ol3

texts:
I am using OL3, Geoserver.
I made a website with Asp.net and add map with OL3. which is get data from geoserver with WFS. my website is located in C:\programfile\geoserver\www\mysite.
Also config IIS that when client (with IP 192.168.0.100) clicl IP of server, it will show the website. (server IP is 192.168.0.1)
the problem is that I receive this error.

XMLHttpRequest cannot load http://localhost:8000/geoserver/cite/ows?service=WFS&version=1.0.0&request=…e&typeName=cite:MyFeature&maxFeatures=50&outputFormat=application/json. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://192.168.0.100' is therefore not allowed access.

how can I config IIS and Geoserver that solve this problem?


texts: Delete duplicate polygon nodes with PostGIS

texts:
I'm trying to delete duplicate nodes from a shapefile with QGIS, but after running the tool to delete it, nothing happens and shows the same errors.

I thought that maybe importing to postgis could be easier. What I want to do is:

Extract points from polygon (this is already done)
Delete the duplicate points ¿with st_equals?
To finish, reshape the points to polygons without the duplicate points.

Anyone know if this methodology would be ok?
I can do the first and second steps, but don't know if number 3 is possible, or if there is a better alternative.

Any clue or suggestion?


texts: Best way to animate radar data in openlayers

texts:
I have a project using OpenLayers which displays radar data on a map.  The radar data itself is just a bunch of polygons each of which is filled with one of 6 possible colors.  The goal is to animate the radar on the map.  Each file contains radar data for a given time and the files are separated by about 5 minutes so my current approach is to loop through all the files and load them one by one into new separate layers.  Once each layer is created it has its visibility set to false and it's added to the map.  I then animate the layers using a timer which turns on the visibility of one layer and turns off the visibility of the preceding layer.  Currently the layers are all vector layers and the data is loaded from KML files although the data files can be pretty much any format that will work best for this project.

The problem with this approach is that once I hit a fairly large time span (around 3 hours or so) of data (equating to approximately 36 layers) the memory consumption gets pretty high (around 250mb).  The final product is supposed to be able to handle up to 18 hours of data in one loop which based on the above number would require over 1GB of memory just for the animation and would likely crash the browser or at least make it very sluggish.

I attempted the same thing using WMS layers for each layer but the redrawing was too slow (the animation changes layers every 100ms) and the memory consumption wasn't much better than the vector version.

I have scoured the net trying to find some example of how to do this but so far am coming up empty handed.  Does anyone out there have any suggestions on how to do this?  I'll take any ideas you have at this point as I've been stuck on this for weeks now with no answer in sight.


texts: Transforming featureCollection back to an actual image


texts: To select the range of values using Raster calculator in QGIS

texts:
I am working on a project where, I need to give ratings for the range of values.
I have the following ratings

index     ratings
0 to 4       4
5 to 8       3
9 to 11      2
12 to 15     1

I am using the syntax as,

(0<= "slope_1@1" <= 4) = 4

The software identifying the syntax as valid but the result showing in null...

please help me from this...


texts: Correct use of buffering by field (distance) in ArcGIS for Desktop?

texts: I am trying to create a buffer in meters around some points using a field with the distance. However I cannot assign a distance type to use and by default I assume it's doing decimal degrees (I am ending up with 'squished' buffers). Is there any way to use meters when buffering by a field?


texts: Counting rocks of different sizes on shoreline using ArcGIS Desktop?

texts:
I have a 5cm resolution geotiff image for a shoreline. I need to count the number of rocks in three classes (e.g., >6 feet, 3-6 feet, and <3 feet). I don't have any point cloud data.

What will be the best gis approach to solve this problem?


texts:
I am looking for ways to remove isolated roads R after filtering out the roads of interest (see attached image). These roads were not isolated in the original .osm file, but became isolated after some filters were applied. They might have been isolated after filtering as they might have been mis-tagged in the original data.

The road data (.osm file) is obtained from Geofabrik, and the desired location selected using a combination of osmosis and using the osmar package in R (see code snippet).

## import osmar
library(osmar)
## Import road.osm file
muc <- get_osm(complete_file(), osmsource_file(file="road.osm"))
## discard paths with value 'footway' as attribute
hw_ids <- find(muc, way(tags(k == "highway" & v!="footway")))
hw_ids <- find_down(muc, way(hw_ids))
hw <- subset(muc, ids = hw_ids)
## plot
plot(hw, col= "black")

Any help using some combination of osmosis and R to remove these isolated roads is preferred. I have seen a solution on gis that suggested the use of ArcGIS, but I do not have access to ArcGIS.


texts: Is there an OSM High Level Dataset (pbf) available?

texts:
I'm currently trying to setup an offline Nominatim geocoder VM and I only need to have access high level data (countries, states and cities).

Is there an OSM high level dataset (.pdf package) available?
A lighter version of Planet.pbf would be perfect for my needs (ref: geofabrik).

Thanks


texts: Explanation of the Thinness ratio formula?

texts:
I am looking for sliver polygons and am using the following formula to identify which polygons have a smaller area-to-circumference ratio (aka Thinness Ratio):

4 * pi * area/(length*length)

That much I understand. But what is not fully clear, is the 4 * Pi bit and why the length has to be squared. Can someone explain this in simple terms?


texts: Error when trying to display a raster from command line! (GRASS GIS)

texts:
I am using GRASS GIS 7.8 on Ubuntu 20.04. Currently working my way through a couple of tutorials. I am trying to display a raster from the command line using the d.rast command on a new monitor that I've opened up before (d.mon=wx0). In the shell there is no Error, however a new window pops up containing the following error:
Rendering failed: Expecting 720x480 image but got 720x456 image.

I am using the spearfish dataset from the GRASS GIS website and the raster I want to display is the spot.image. Any clew what might be the problem?


texts: vb.net arcmap add-in loop through a dbf file get 'System.NullReferenceException' error

texts:
I am using vb.net arcmap add-in.

I encounter a System.NullReferenceException error and cause arcmap crash.
The workflow in my vb.net code is:

(1) create a ShapefileWorkspaceFactory

(2) use the workspace to open a DBF file as ITable

(3) loop through each record and get related value

The code I am using:

    Dim pWSF As ShapefileWorkspaceFactory = New ESRI.ArcGIS.DataSourcesFile.ShapefileWorkspaceFactory()
    Dim pWS As ESRI.ArcGIS.Geodatabase.IWorkspace = pWSF.OpenFromFile("D:/dev/local", 0)
    Dim pFws As ESRI.ArcGIS.Geodatabase.IFeatureWorkspace = CType(pWS, ESRI.ArcGIS.Geodatabase.IFeatureWorkspace)
    Dim pTable As ESRI.ArcGIS.Geodatabase.ITable = pFws.OpenTable("Segment_line_intersect")
    Dim pCur As ESRI.ArcGIS.Geodatabase.ICursor = pTable.Search(Nothing, False)
    Dim queryFilter As IQueryFilter = New QueryFilterClass()
    Dim pFieldIndex1 As Integer = pTable.Fields.FindField("Segment_ID")
    Dim pFieldIndex2 As Integer = pTable.Fields.FindField("DrainName")

    Dim pRow As ESRI.ArcGIS.Geodatabase.IRow = pCur.NextRow
    Dim SegmentID As String
    Dim DrainName As String

    Do Until pRow Is Nothing
        pRow = pCur.NextRow

        SegmentID = pRow.Value(pFieldIndex1)
        DrainName = pRow.Value(pFieldIndex2)

        MsgBox("SegmentID = " + SegmentID)
        MsgBox("DrainName = " + DrainName)

    Loop

All the values of SegmentID and DrainName can be retrieved correctly, but once all the values retrieved, the map crash. I use debug mode to track the error and get 'System.NullReferenceException' error.

I believe there is something wrong with my code, but I cannot tell what cause the issue.

I am using arcmap 10.2.2.


texts: What is appropriate coordinate transformer for Lambert Azimuthal Equal Area to Google maps in Geotools?

texts:
I need to extract the coordinates from the GeoTiff file (that was provided by European Soil Agency for Spain) and display them in Google Maps.

For that, I wrote simple Java code

    public static void readTiff1(final String pathToTiff) {
        var file = new File(URLDecoder.decode(TiffParser.class.getResource(pathToTiff).getFile(), StandardCharsets.UTF_8));
        var format = GridFormatFinder.findFormat(file);
        var reader = format.getReader(file);
        var cov = reader.read(null);
        var gridRange2D = cov.getGridGeometry().getGridRange();
        var googleMapsTransformer = CRS.decode("EPSG:3857");
        var target = cov.getCoordinateReferenceSystem();
        var gf = new GeometryFactory();
        var targetToWgs = CRS.findMathTransform(target, googleMapsTransformer, true);
        int count = 0;
        for (int j = gridRange2D.getLow(1); j < gridRange2D.getHigh(1); j++) {
            for (int i = gridRange2D.getLow(0); i < gridRange2D.getHigh(0); i++) {
                if (count++ > 10)
                    System.exit(1);
                var coord = new GridCoordinates2D(i, j);
                var p = cov.getGridGeometry().gridToWorld(coord);
                var point = gf.createPoint(new Coordinate(p.getOrdinate(0), p.getOrdinate(1)));
                var wgsP = JTS.transform(point, targetToWgs);
                System.out.format("(%d %d) -> POINT(%.2f %.2f) -> POINT(%.2f %.2f)%n", i, j, point.getCoordinate().x,
                        point.getCoordinate().y, wgsP.getCentroid().getCoordinate().x, wgsP.getCentroid().getCoordinate().y);
            }
        }
    }

The output of this program is the following

PROJCS["PCS Name = GRS_1980_IUGG_1980_Lambert_Azimuthal_Equal_Area",
  GEOGCS["GCS Name = GCS_GRS_1980_IUGG_1980",
    DATUM["Datum = unknown",
      SPHEROID["Ellipsoid = GRS80", 6378137.0, 298.257222101]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH]],
  PROJECTION["Lambert_Azimuthal_Equal_Area"],
  PARAMETER["latitude_of_center", 51.99999999999999],
  PARAMETER["longitude_of_center", 10.0],
  PARAMETER["false_easting", 4321000.0],
  PARAMETER["false_northing", 3210000.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH]]

And the list of the coordinates

indexes  source                          target
(0 0) -> POINT(2635950.00 5421950.00) -> POINT(-3503460.95 10202813.83)
(1 0) -> POINT(2636450.00 5421950.00) -> POINT(-3502415.88 10203547.81)
(2 0) -> POINT(2636950.00 5421950.00) -> POINT(-3501370.61 10204281.72)
(3 0) -> POINT(2637450.00 5421950.00) -> POINT(-3500325.13 10205015.57)
(4 0) -> POINT(2637950.00 5421950.00) -> POINT(-3499279.46 10205749.36)
(5 0) -> POINT(2638450.00 5421950.00) -> POINT(-3498233.59 10206483.08)
(6 0) -> POINT(2638950.00 5421950.00) -> POINT(-3497187.51 10207216.75)
(7 0) -> POINT(2639450.00 5421950.00) -> POINT(-3496141.24 10207950.34)
(8 0) -> POINT(2639950.00 5421950.00) -> POINT(-3495094.76 10208683.88)
(9 0) -> POINT(2640450.00 5421950.00) -> POINT(-3494048.09 10209417.35)
(10 0) -> POINT(2640950.00 5421950.00) -> POINT(-3493001.21 10210150.76)

What is the appropriate coordinate transformer for the google maps api?


texts: Export one point with shared origin value in attribute table

texts:
I am trying to find a way to export unique points from a pointlayer in QGIS. Using the tool ‘Extract vertices’, I generated corner points from polygons in the input layer (first image). I now only need one point per polygon to remain in the pointlayer. The ‘Extract vertices’ tool conveniently generates a column in the attribute table containing the vertex_part_ring value, portraying from what vertex the point is originated from (second image). I am looking for a way to (automatically) export one point from a shared vertex_part_ring series.
Below is the polygon the points are originated from. (Sample project)

Below is a screen capture of the attribute table as a result of the extract vertices tool. The column ‘vertex_part_ring’  gives the origin of the point. I need to keep only one point from each origin series.


texts: Is there a way to remove this thick gray border from my raster symbology in the legend?

texts:
I have this problem with different raster layers, in this instance it's land use data. In my layout, and in a rendered PDF of my map, the legend entry for the raster layer has an ugly, thick gray border around my raster color patches. Here's how they look when rendered to PDF:

Sometimes it's even worse, and ArcGIS appears to randomly choose a thick border for some colors and a thinner, lighter one for others, making my legend very ugly. Here's another example in a finished PDF:

Here's what the same legend looks like in the ArcGIS layout. The border width looks the same, nothing like the PDF. The color differences are not as pronounced as in the PDF but you can see they're there:

Is there any way to edit or remove those gray borders? I've searched a bit online and the most recent thread I've found is this one from the ArcGIS forums, with a couple workarounds that I may resort to, but I'm hoping for a simpler solution.

FYI, my rasters' symbology is set to use Unique Values, but using Colormap gives the same results.


texts: Relative Hyperlinks in QGIS without end user's software dependency

texts:
I am trying to work out how can I add a relative path hyperlink to a shapefile in QGIS. I have found several great tips on various forums but only one has so far potential to be implemented. The overall goal is for the end user to be able to access a hyperlinked pdf/png/jpg etc without me having to set up any software dependencies under layer properties Actions tab in QGIS. I assume that the end user will have the software required to open a document but I would like the end user's computer to choose a default program to open the files. One of the possibilities would be to use a 'Python' action which I located on one of the forums which after several modifications did not work for me: http://lists.osgeo.org/pipermail/qgis-user/2011-May/012117.html.

Here is the screenshot of my layer's action properties:

In the attribute table, the column with hyperlink is called 'PDF and the relative hyperlink is: /doc/document.pdf

My folder structure is:

 - /HyperlinksTest
 -- /data
 --- Renewable_Energy_Zones.shp
 -- /doc
 --- document.pdf
 HyperlinksTesting.qgis

When I identify the feature, I get the following Python error:

Am I using wrong variables in the Python code?
Is there any other way of using relative paths in hyperlinking data to external documents? I noticed there is eVis extension but as far as I see it relies on specyfying software to open a document with, which I'd like to avoid.

Many thanks,
Magda


texts:
I'm using ArcGIS Server 10.0. I have a python geoprocessing that takes a lot of time. The geoprocessing calculates geometries using points, polygons, buffers, and cursors.

I'm also using sql queries with the SdeExecuteSQL function. This function seems to take time too. Is it possible to optimize this?

Do you have any hint of how to optimize the geoprocessing?

Would it be shorter in time if I used the Java API?


texts: Calculating Runoff / Streams using DEM in ArcGIS

texts:
I am just curious about the workflow that many of you use for runoff/stream analysis using a DEM and create a "water network" from the results in ArcGIS 10.1? I know the average rainfall for different seasons throughout the year, but that's really all I have to work with for my area.

I mainly want to investigate where water will "end up" given a fairly large rain event. I am guessing you would use ArcHydro for this? Is there a different software suite or suite of tools that you would recommend for this?

Thanks!


texts: Consuming ArcGIS MapServer as FeatureServer in ArcGIS Desktop?

texts:
I work on ArcGIS 10.3, and I try to find a way to convert MapServer into FeatureServer.

The MapServer has no attribute table, and there no way to change the symbology.
The FeatureServer has attribute table and after I add it to the map, i can export it into feature layer and it acts as a regular layer.

I found the layer here:

https://gwa.israntique.org.il/mobidic/rest/services

Is there any way to do it?


texts:
I work for a small company that designs alignments using ArcGIS.  (To give you an idea of company size, we have about 50 employees; I am the only full-time GIS person, although we have 4 ArcMap licenses.)

In the past, the alignments we have dealt with have been fairly straightforward (5 miles or less) so we have handled the numerous version of the alignments simply by saving dated shapefiles for each version of the alignment.

Now we are taking on a project that will be significantly longer and more complex (30+ miles) so I am looking for ways to track the version and changes to our alignment.

I have read some about ArcMap's versioning capabilities but this does not seem like exactly what we are looking for.  I do not need the ability for multiple users to edit the database; I will be the only user making changes to this database.  I do not have any kind of ArcSDE geodatabase set up.

What I really need is:

A way to keep track of each version of the alignment, so I can see what it looked like at any particular moment in time, and
A way to track/tag each specific change from the previous version (such as the alignment was moved to the other side of a road at the certain location due to overhead utilities)

Can anyone suggest an ArcMap tool or other system or procedure for efficiently keeping track of these kind of changes?

Note: We only have the Basic ArcMap License.  Geodatabase Archiving seems like a great tool for this use, but it requires a Standard or Advanced level license.  Any suggestions for tools that work with an ArcMap Basic license?  Upgrading would be an option to look into, but I want to explore other options first since Standard may be out of our price range.