texts: Using true north bearing in label expression

texts:
In QGIS there is an option to rotate the north arrow as grid north or true north. For my case I need the true north angle. The rotation works just fine but now I also wanted to add the image rotation of said north arrow as a text label. Now I wondered if I can query the value which is actually displayed inside the item properties of the north arrow.

Up to now I tried grabbing it using item_variables('northarrow') with no success. Furthermore I tried using a python function in the label expression utilizing bearingTrueNorth() which looks as following:

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

@qgsfunction(args='auto', group='Custom')

def truenorth(value1,feature,parent):

 crs=QgsCoordinateReferenceSystem.fromEpsgId(4326)
 transformContext=QgsCoordinateTransformContext()

 bearing=bearingTrueNorth(crs,transformContext, QgsPointXY(value1))

return bearing

I tried various iterations but to no avail. My value1 is as QgsGeometry grabbed by using centroid(@atlas_geometry).

Is there any known solution to this?


texts:
I've worked a little with basic exporting of CAD .DXF files from our software using GDAL.  We want to do more sophisticated things, like importing and creating georeferenced files with multiple layers, and it is clear that GDAL's functionality around CAD formats is rather limited.  Other questions like this one address the lack of .DWG format support for various reasons (explained in more detail at GDAL's website).  There is good information about open source CAD software in other questions like this.

Does anyone know of good open-source software libraries other than GDAL that have more sophisticated CAD file format reading and creation abilities we could use in our software to create .DXF or .DWG files?

Should such a library not exist, an alternative workflow like using file conversion sounds possible, like using FME to convert from different formats to CAD as discussed in this post.


texts: Kriging using multiple spatial input variables and one spatial output (response) variable

texts:
I have a data table with the following columns in a CSV file.

X    Y     Zinc      L    F     C
---  ---   -----     --   ---  ----

Of course, I can import the CSV file as a data frame

a <- read.csv('givendata.csv', header = T)

And convert this data into spatial data by

coordinates(a) = ~X + Y

Now, I wish to run kriging by using L, F and C as spatial input variables and Zinc as the output variable. How do I achieve this?

Edit: I am new to Kriging.

I am adding some more explanation to the question, based on a comment that I do not have an idea about linear modelling.

If I do the following

A1 <- read.csv('givendata.csv', header = T)
mymodel <- lm(Zinc ~ L + F + C, data = A1)

it gives the output, with all the regression coefficients and the constant. But, these observations have spatial inputs (see the meuse dataset), which is why it is referred to as a spatial variable (this is GIS.SE!).

Now suppose I were to take the following steps:

x.range <- (range(A1@coords[,1]))
y.range <- (range(A1@coords[,2]))

A1.grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=30),
                    y=seq(from=y.range[1], to=y.range[2], by=30))

coordinates(A1.grd) <- ~x+y
gridded(A1.grd) <- TRUE

library(automap)
kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd)

I get the following error:

Error in eval(expr, envir, enclos) : object 'L' not found

which is funny because the lm was working and the Kriging does not work!


texts: MapFish print module with cql_filter params

texts:
Does anyone know how to pass cql_filter params into MapFish's Print module?
I already tried this, but it does not work. It generates the PDF but I get this error "The byte array is not a recognized imageformat.
":

{
  "units" : "degrees",
  "srs" : "EPSG:4326",
  "layout" : "A4 portrait",
  "dpi" : 75,
  "mapTitle" : "Printing Demo",
  "comment" : "This is a simple map printed from GeoExt.",
  "layers" : [{
      "baseURL" : "http://zzz.com:8080/geoserver/wms",
      "opacity" : 1,
      "singleTile" : true,
      "type" : "WMS",
      "layers" : ["zzz:yyy"],
      "customParams" : {
        "CQL_FILTER" : "dico in('1106')"
      },
      "format" : "image/jpeg",
      "styles" : [""]
    }
  ],
  "pages" : [{
      "center" : [-9.157682, 38.741626],
      "scale" : 4000000,
      "rotation" : 0
    }
  ],
  "legends" : [{
      "name" : "Tasmania",
      "classes" : [{
          "name" : "",
          "icons" : ["http://zzz.com:8080/geoserver/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetLegendGraphic&EXCEPTIONS=application/vnd.ogc.se_xml&LAYER=zzz:yyy&FORMAT=image/gif&SCALE=4000000"]
        }
      ]
    }
  ]
}


texts: I am trying to make a DEM using topo to raster tool. When it is done a message pops up that some drawings Cannot acquire a lock, and the DEM is not visible. I have made a DEM before with the same data and it worked.


texts: Style layer on server side or client side?

texts: I have GeoServer installed where I am planning to place our company layers. Since we are using a web front end (Leaflet) to get data onclick (GetFeatureInfo) from some WMS layers should I style the GetFeatureInfo response on server side or client side?


texts: How can I handle ArcMap events of CREATE, UPDATE or DELETE?

texts: As a developer I can use ArcPy (Python) or ArcObjects (.NET) to develop some extensions or tools for ArcGIS platform. But I need to handle some events. When a user connected to an Enterprise geodatabase and create an object in database table, can I handle the events of CREATE, UPDATE or DELETE? Is this possible?


texts: Reverve polyline direction in mapinfo 15

texts:
How can I reverse the line direction in MI 15?

I have tools for older versions but can't find the same for the current version - http://www.mapinfotools.com/index.php?option=com_sobi2&sobi2Task=sobi2Details&catid=10&sobi2Id=329&Itemid=53 and http://www.geoas.net/Help/1200/GeoAS_MapCAD_Pro_en/index.html?reverse_line_direction.htm


texts:
I need to use arcpy.Select_analysis() in my arcpy script, but I have problem with using variables. When I use hard coded values in where_clause, it is working.

where_clause = '"Area" > 20'
arcpy.Select_analysis("input", "output", where_clause)

When I try to use variables in where_clause because user of toolbox will enter name of field and value of minimimum area, there is always problem with invalid expression. I tried lots of types of expressions, but nothing is working.

field_area = "Area"
minimum_area = 20
where_clause = 'field_area > minimum_area'
arcpy.Select_analysis("input", "output", where_clause)

Error:

  Runtime error  Traceback (most recent call last):   File "",
  line 1, in    File "c:\program files
  (x86)\arcgis\desktop10.2\arcpy\arcpy\analysis.py", line 98, in Select
  raise e ExecuteError: ERROR 000358: Invalid expression field_area >
  minimum_area Failed to execute (Select).

I work in ArcMap 10.2.


texts: How to show label conditionally on the value of an attribute?

texts:
In the MapServer Mapfile, I have this:
LABEL
    TEXT "Q=[Q] m³/h"
END

The problem is that for those features having an empty (null) Q, the result is Q= m³/h. Instead, I'd like it to not show anything. How can I achieve that?
Can I also conditionally exclude parts of the string? E.g.
LABEL
    TEXT "Q=[Q] m³/h D=[Diameter] mm"
END

Here the Q=[Q] m³/h part should show only if Q is not null, and the rest only if Diameter is not null.


texts:
I am trying to convert a GeoTIFF to a KMZ. Here's my code:

from osgeo import gdal
import zipfile

srcf = r"C:\raster_to_kml\OrthoImage.tif"
vfn = "/vsimem/OrthoImage.jpg"
drv = gdal.GetDriverByName('JPEG')
src = gdal.Open(srcf)
tgt = drv.CreateCopy(vfn,src)
rasterLyr = QgsRasterLayer(srcf,"OrthoImage")
e = rasterLyr.extent()

kml = """<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
    <Document>
        <name>Valleyfield</name>
        <GroundOverlay>
            <name>Image</name>
            <drawOrder>30</drawOrder>
            <Icon>
             <href>OrthoImage.jpg</href>
            </Icon>
            <LatLonBox>
              <north>%s</north>
              <south>%s</south>
              <east>%s</east>
              <west>%s<west>
            </LatLonBox>
        </GroundOverlay>
    </Document>
</kml>""" %(e.yMaximum(),e.yMinimum(), e.xMaximum(), e.xMinimum())

vsifile = gdal.VSIFOpenL(vfn,'r')
gdal.VSIFSeekL(vsifile, 0, 2)
vsileng = gdal.VSIFTellL (vsifile)
gdal.VSIFSeekL(vsifile, 0, 0)
z = zipfile.ZipFile(r"C:\raster_to_kml\OrthoImage.kmz", "w", zipfile.ZIP_DEFLATED)
z.writestr("doc.kml", kml)
z.writestr("OrthoImage.jpg", gdal.VSIFReadL(1, vsileng, vsifile))
z.close()

I use QGIS 3.0
Here's the error message i get:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 38, in <module>
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\zipfile.py", line 1760, in writestr
    zinfo.file_size = len(data)            # Uncompressed size
TypeError: object of type 'NoneType' has no len()


texts: ArcGIS Android Reverse geocoding, Cannot perform query. Invalid query parameters

texts:
I'm trying to write reverse geocoding. And i copied code given here... Here is my code sample..

      public class MagProActivity extends Activity {

            private MapView mMapView;
            private GraphicsLayer locationLayer;

            /** Called when the activity is first created. */
            @Override
            public void onCreate(Bundle savedInstanceState) {
                super.onCreate(savedInstanceState);
                requestWindowFeature(Window.FEATURE_NO_TITLE);
                setContentView(R.layout.main);

                mMapView = (MapView) findViewById(R.id.map);
                mMapView.addLayer(new ArcGISDynamicMapServiceLayer("http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"));

mMapView.setOnLongPressListener(new OnLongPressListener() {

            private static final long serialVersionUID = 1L;

            @Override
            public void onLongPress(float x, float y) {
                Point mapPoint = mMapView.toMapPoint(x, y);

                new ReverseGeocoding(MagProActivity.this, mMapView).execute(mapPoint);

            }
        });

      }
    }

ReverseGeocoding class

import java.util.Map;

import android.app.Activity;
import android.os.AsyncTask;
import android.widget.Toast;

import com.esri.android.map.MapView;
import com.esri.core.geometry.Point;
import com.esri.core.geometry.SpatialReference;
import com.esri.core.tasks.ags.geocode.Locator;
import com.esri.core.tasks.ags.geocode.LocatorReverseGeocodeResult;

public class ReverseGeocoding extends AsyncTask<Point, Void, LocatorReverseGeocodeResult>{

    private Activity curActivity;
    private MapView mapView;

    public ReverseGeocoding(Activity _curActivity, MapView _mapView){
        this.curActivity = _curActivity;
        this.mapView = _mapView;
    }

    @Override
    protected LocatorReverseGeocodeResult doInBackground(Point... params) {
        // create results object and set to null
        LocatorReverseGeocodeResult result = null;
        // set the geocode service
        Locator locator = new Locator();
        try {

            // Attempt to reverse geocode the point.
            // Our input and output spatial reference will be the same as the map.
            SpatialReference mapRef = mapView.getSpatialReference();
            result = locator.reverseGeocode(params[0], 50.0, mapRef, mapRef);
        } catch (Exception e) {
            e.printStackTrace();
        }
        // return the resulting point(s)
        return result;
    }

    protected void onPostExecute(LocatorReverseGeocodeResult result) {

        String resultAddress = "";

        // Construct a nicely formatted address from the results
        StringBuilder address = new StringBuilder();
        if (result != null && result.getAddressFields() != null) {
            Map addressFields = result.getAddressFields();
            address.append(String.format("%s\n%s, %s %s",
                    addressFields.get("Street"), addressFields.get("City"),
                    addressFields.get("State"), addressFields.get("ZIP")));

            // Show the results of the reverse geocoding in a toast.
            resultAddress = address.toString();
            Toast.makeText(curActivity, resultAddress, Toast.LENGTH_LONG).show();
        }
    }

}

When i long click on mapview this error will come com.esri.core.io.EsriServiceException: Cannot perform query. Invalid query parameters and toast will show all parameters null except city.. I think i'm doing some mess with Tiled service on Dynamic service.. Please help what i did wrong..

Thanks in Advance


texts: I wish to use rounded values as labels using QGIS 3.4.4. My field is type "double" and the the data source is a Spatialite table. I notice two odd behaviors linked to this problem. First, the attribute values have false precision. For example, in the source file, a value is 16.300000, but it appears as 16.300000000000001 in QGIS. Second, when I attempt to use the round function under labels, the value can either be truncated to zero decimal places using round or the precision can't be reduced at all, no matter what precision I supply to round. How can I resolve this problem?


texts: Alternative to Snap (Editing) in QGIS/MapInfo?

texts: Is there any alternative to Snap (Editing) in either QGIS or MapInfo (user-made tools), or another way to do this in ArcGIS without an advanced license?


texts: Setup GeoNetwork and Geoserver withTomcat

texts:
I need to set up GeoNetwork on the same server as GeoServer (Windows server 2003, with IIS6) .

Is it best to use multiple instances of Tomcat or put both web archives in the same one? They both need to listen on the same port 8080.

Thanks


texts:
I am using Landsat Surface Reflectance Level-2 images. I have to convert the DNs to reflectance values. When I apply the formula provided in the Landsat8 User Handbook :

  reflectance = Reflectance_MULT_Band*DN  + Reflectance_ADD_Band

in the metadata file, Reflectance_MULT_Band = 2.0000E-05 and Reflectance_Add_Band= 0.1
The problem is I am getting negative values even after  correction for
the solar elevation angle.


texts: Summing values in one field using PyQGIS

texts:
In QGIS 2.18 there is a shapefile with three fields accordingly

id
Type
Value

1
a
5

2
b
NULL

3
c
1

4
d
6

5
e
NULL

I want to get the total sum of all values in the field "Value". How can I achieve that in the Python Console?
The field "Value" maintains two types of data int and NULL.
I was trying with a while-loop but it does not work
layer = qgis.utils.iface.activeLayer()

if not layer.isValid():
    raise Exception('Layer is not valid')

features = layer.getFeatures()

for feat in features:
    attr = feat.attributes()[2]

    if attr == NULL:
        continue

    else:
        total = 0
        index = 0
        while index < layer.featureCount():
            total = total + int(attr)
            index = index + 1
return(total)

How to achieve one number that gives the total sum of a whole field exclusive NULL-values?

After executing @ahmadhanb's solution I received a correct number but
somehow with a strange output type, when print (type(sum(total))) it gives me
<type 'NoneType'>. What could be a problem?

References:

Find Neighbor Polygons in a Layer
Summing values from different intersecting/overlapping points in QGIS
Sum new values to existing attribute values in PyQGIS?
SQL query in Python console to sum column based on values of other two columns


texts: In OpenLayers how can I enable users to cancel the completion of a polygon, and carry on drawing it?

texts:
I am using the OpenLayers v4.6.5 Draw interaction.

drawInteraction = new ol.interaction.Draw({...});
map.addInteraction(drawInteraction);

I want users to be able to add properties to the features they draw so use the drawend event to display a form.

drawInteraction.on('drawend', function(e) {...}).

This all works fine.

However, an issue that arises is that, when drawing polygons, users sometimes complete polygons without intending to (for example by clicking within the snap distance of the start point or clicking too fast so triggering a double click event). On simple polygons they can easily cancel and start again, but many of our polygons are complex.

I wondered if there is a way, once the form is displayed, of undoing their last click and cancelling the completion of the geometry? I imagine them clicking an "Undo last click" button on the form and just carrying on with drawing their polygon. I have explored removeLastPoint (but this seems to have no effect once the drawing has finished) and finishCondition (but cannot see how to use it achieve what I want).

Users are able to modify their incorrect polygons to get them to be what they want, so it is not a huge issue. I understand I could make it harder for them to accidentally complete their polygons, by changing snap distances etc.

But an undo functionality would make it more intuitive for them. How can it be done?


texts: How do I reorder the Value field in a raster? There is no up or down button in it.


texts: Splitting lines/creating vertex at intersection using PyQGIS?

texts:
I am looking for an easy way or function to split 2 previously selected line features (generate vertex) at intersecting points (or touching points). See the following picture:

I am developing a function/plugin where the user is supposed to manually select 2 line features which intersect and wants them to be splitted at intersection (touching) point.

So something like this:

layer = iface.activeLayer()
feats =  layer.selectedFeatures()

geom1 = feats[0].geometry()
geom2 = feats[1].geometry()

# I need a function like this
geom1.splitAtIntersection(geom2)

The line features/layer are loaded from a PostGIS database, so it is important to edit the existing layer and not create a new layer or shapefile. It would be nice if there is a pure pyqgis solution which I overlooked to solve this problem.

I found the intersects() and touches() functions of the QgsGeometry class, but they dont return the intersection points. There has to be a native pyqgis function to retrieve the intersection point or to split lines at a given point, right?

I know there is the split lines with lines tool, but I only want to split manually selected features/lines and not all.

I am aware that there are several similar questions asked already, but none of them seem to solve my particular issue.

For example:

Intersecting lines get crossings - no pure PyQGIS solution (using shapely, fiona) and generating a new shapefile.
Create vertex where two lines meet - generating a new layer which I do not want (rather edit existing)


texts: Displaying rtmp live stream video in Cesium?

texts:
I would like to display an rtmp stream video in Cesium. I created an entity and used the <video> element as its material. When I use video stored locally (in .mp4 format) it worked fine, but when I switched to an rtmp source and play it via videojs, it can play in <video> element, but it can't play on the entity. The rtmp uses flash to play, while ordinary video does not require flash, I wonder if its because of this.

Here's the <video> element

<video id="rtmpStream" class="video-js" controls preload="auto" width="640" height="300" data-setup="{}" autoplay>
        <source src="rtmp://202.69.69.180:443/webcast/bshdlive-pc" type="rtmp/flv">
    </video>


texts: How to perform regression Kriging with R studio using gstat?

texts:
I'm struggling with regression kriging in R. I want to interpolate the temperature of a watershed under consideration of elevation, but get always the same error when I want to perform the last step, which is the kriging.
I have two data frames, one of my 15 gauging stations with coordinates, elevation and the meteorological data for around 30 years on a daily basis, each day is one column (X1....X11642). The other data frame contains information about my watershed as point data with coordinates and elevation (as grid code in table)(6 columns, 315'226 rows).
My code looks like this:
library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
library(gstat)

d_t<-data_temperature
coordinates(d_t)= ~LONG+LAT

jh_sp<- jh_coord
coordinates(jh_sp)=~X+Y

ext_jh<-extent(jh_sp)

jh_ras<- raster(ext_jh)

jh_raster<-rasterize(jh_sp, jh_ras,jh_sp$GRID_CODE, fun='last', background=NA,
          mask=FALSE, update=FALSE, updateValue='all', na.rm=TRUE)

proj4string(jh_raster)<-CRS("+init=epsg:4326")

jh_grid <- rasterToPoints(jh_raster, spatial = TRUE)

gridded(jh_grid)<- TRUE

proj4string(d_t) <- CRS("+init=epsg:4326")

stat_prj_t <- spTransform(d_t, CRSobj = CRS("+init=epsg:4326"))

lzn.vgm = variogram(log(X1)~1, stat_prj_t)

lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 4000, 1))

#plot(lzn.vgm, lzn.fit)

lzn.kriged = krige(stat_prj_t$X1~stat_prj_t$ELEVATION , newdata=stat_prj_t, locations=jh_grid, model = lzn.fit)

#gjh <- gstat(id = "stat_prj_t", formula = stat_prj_t$X1~stat_prj_t$ELEVATION,
                #data = stat_prj_t, maxdist=100, nmin=10, force=TRUE, model=lzn.fit)

#predgjh<- predict(gjh, jh_grid)

#spplot(predgjh)

When I perform:
lzn.kriged = krige(stat_prj_t$X1~stat_prj_t$ELEVATION , newdata=stat_prj_t, locations=jh_grid, model = lzn.fit)

I get following error:
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  :
  dimensions do not match: locations 154 and data 15

When I perform:
gjh <- gstat(id = "stat_prj_t", formula = stat_prj_t$X1~stat_prj_t$ELEVATION,
                data = stat_prj_t, maxdist=100, nmin=10, force=TRUE, model=lzn.fit)
predgjh<- predict(gjh, jh_grid)

I get this error:
Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  :
  NROW(locs) != NROW(X): this should not occur

I tried to figure out how to solve one of these problems but wasn't successful, yet.
I already searched for answers on google and here in the forum, but non of suggestions worked for me.
How can I perform regression kriging without problems in R studio?


texts: How to add data to a newly created folder using arcpy?

texts:
I'm writing a script to create a new folder and then copy a feature class to a shapefile and add the shapefile in the new folder. The folder gets created but the script is failing on the feature class to shapefile conversion line with an error of:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Output Folder: Dataset Zoning_ does not exist or is not supported
Failed to execute (FeatureClassToShapefile).

My code so far is:

import arcpy, os, datetime
from arcpy import env

# Set Workspace
env.workspace = "C:/GIS_Work"

inputFile = arcpy.GetParameterAsText(0)
outputFile = arcpy.GetParameterAsText(1)

# Create folder with today's date
#dt = str(datetime.datetime.now())

# Set local variables
out_folder_path = "C:/GIS_Work/PropertyAppraiser"
out_name = "Zoning_"

# Execute CreateFolder
arcpy.CreateFolder_management(out_folder_path, out_name)

# Execute Export Feature Class to Shapefile
arcpy.FeatureClassToShapefile_conversion(inputFile, out_name)

Is there something I need to do to the folder to make it able to receive files? I couldn't see anything in the ArcGIS Help that would suggest this, but since the folder is getting created the portion of the error that says "does not exist" wouldn't apply, I would think.


texts: Get area and perimeter of 3D plane using JTS library

texts:
is there any Java API in Geotools JTS lib for calculating area of 3D plane?
Same about Line, where points have three dimensions.


texts: ArcGIS raster calculator: ERROR 000539: Error running expression: rcexec() : No module named collections

texts:
My raster is in the GRID format, source type is continuous and pixel type is floating point (32bit).

If i use the raster calculator ("dem" *1000 or "dem >= 2m" (to highlight certain heights) it gives me the following error:

  ERROR 000539: Error running expression: rcexec() : No module named collections.


texts: Mars DTM has all negative elevations in QGIS. Can't build STL file for 3D print

texts:
I downloaded Mars HiRise DTM .IMG files and they load fine and generate contours fine. When I try to create an STL using DEMto3D to print it on the 3D printer I can't export it because I cannot put any value in the height field... I get an error if I put in any positive value that that says:

Height of the base must be lower than DEM highest point.

Well the highest point is a NEGATIVE altitude. The DTM goes from -3190 up to -2050... As soon as you put a negative sign in the height field it generates a "ValueError:could not convert string to float"

There is no good documentation at all on the DEMto3D plugin so I have to ask here if I can somehow fool the data to get an STL out of this?


texts:
'NV' cells in a raster layer appear to be exported to a .png as transparent (i.e. white), despite being recoloured in GRASS. Is there any way to force the cells to keep their colour?

Edit: A workaround is to open the .png in paint and fill in the null values manually


texts: Is it possible to adjust the length of a GeoJSON LineString?

texts:
Let's say I have the following GeoJSON:

{
    "type": "FeatureCollection",
    "features": [{
            "type": "Feature",
            "properties": {},
            "geometry": {
                "type": "Point",
                "coordinates": [-88.0344009399414, 43.11040115356445]
            }
        },
        {
            "type": "Feature",
            "properties": {},
            "geometry": {
                "type": "Point",
                "coordinates": [-89.3375015258789, 43.13990020751953]
            }
        },
        {
            "type": "Feature",
            "properties": {
                "feature_type": "line"
            },
            "geometry": {
                "type": "LineString",
                "coordinates": [
                    [-88.0344009399414, 43.11040115356445],
                    [-89.3375015258789, 43.13990020751953]
                ]
            }
        }
    ]
}

The above is a line between two points. I want to shorten that line by some x amount so it doesn't go all the way to the points at the beginning and end. Is there a way to do this without changing the start and end points for the line? To say, for example shorten the line by 50px or something like that?

I know you can do this when creating a polyline in leaflet, but can't find any settings on a GeoJSON LineString that would allow me to do the same.


texts:
I am completely new to GeoData so this might be somewhat of a silly problem, but I am having a problem with the angle GeoTIFF images tend to come with. I am collecting satellite images for a deep learning segmentation project. For that purpose I need square or rectangular images. But when I import an image with GDAL (via Python) it is an angled image, since these GeoTIFF images contain a projection (in my case an unnamed projection, more info below). When I plot these images I get a rotated images, with the sides filled with black to form a square.

These rotated images are useful in many geo-related cases, but not for deep learning. Every image needs to be a square/rectangle. Every pixel needs to come from the image itself (not a black fill). I want to remove the angle in the image: in essence remove the projection. I have been searching and trying for days to simple remove the projection, but I keep failing. Does anyone here know how I can "unangle" these images, either by removing the projection (if that is even possible) or rotating the image? Or is there an even easier method?

Some information of images I'd like to process (un-angle) is:

Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Amersfoort",
        DATUM["Amersfoort",
            ELLIPSOID["Bessel 1841",6377397.155,299.152812800003,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4289]],

Corner Coordinates:
Upper Left  (  185042.500,  364840.000) (  5d49' 5.55"E, 51d16'22.64"N)
Lower Left  (  185042.500,  348993.000) (  5d49' 0.74"E, 51d 7'49.79"N)
Upper Right (  200000.000,  364840.000) (  6d 1'57.25"E, 51d16'19.08"N)
Lower Right (  200000.000,  348993.000) (  6d 1'50.05"E, 51d 7'46.24"N)
Center      (  192521.250,  356916.500) (  5d55'28.40"E, 51d12' 4.62"N)

This problem does not occur when I import raw satellite imagery, and therefore I am looking in that direction for potential answers. These raw images have a SpatialReference of None, a GeoTransform of (0.0, 1.0, 0.0, 0.0, 0.0, 1.0), and the Corner Coordinates are specified by simple pixel locations (6395 is the size of Y whereas 5720 is the size of X):

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 6395.0)
Upper Right ( 5720.0,    0.0)
Lower Right ( 5720.0, 6395.0)
Center      ( 2860.0, 3197.5)


texts: How to Buffer with Dissolve ALL on very large polyline dataset?

texts:
So I have a fairly large polyline network covering a decent chunk of North America, containing just shy of 900,000 table records. I want to do some very simple "GIS 101" kind of processing to it -- I want to buffer everything with the "Dissolve" parameter set to "ALL".

Well, unsurprisingly, this is way more than ArcGIS is ready to take on, even w/ 6 GB RAM & with everything in a geodatabase. I know the "Dissolve" operation is more prone to running out of memory in such cases than "Buffer", so I tried to run "Buffer" without "Dissolve = ALL", but even that will not generate a usable output (the file exists but won't draw or reveal its attribute table). And, in any case, I am sure the subsequent "Dissolve" would hang.

Here's what I've tried:

Breaking up the network into smaller pieces. Well, this seems like a non-starter, because I need to make sure I don't cut any line segments apart in doing so. Of course selecting polyline segments sequentially from the attribute table
will just select random pieces of polylines all over the place, shredding the integrity of the network.  I'm not sure how to break it up programmatically such that no connected pieces are broken apart.
Running the buffer & dissolve in QGIS, based on some hearsay that it was better at handling large datasets. Like ArcGIS, it also just hung there forever.
Importing everything into R with the "maptools" package and seeing if I could do it there. The data is so huge I couldn't even read it in within a reasonable time span.

There's a whole lot of talk out there about this kind of situation, but it always ultimately ends with "you need to find a way to break up your large dataset." I'm not really sure how this can be done cleanly with a huge network of stream polylines, barring some crazy Python script that can tell, for each polyline segment, which ones it's connected to, and select those only, until the entire spatially discrete sub-network is selected, write that to a single feature, & then move on to the next such sub-network.

This seems absurdly complex as a workaround for very simple processing on a not-even-all-that-big dataset (though I'll try it if I have to), but can anyone tell me how I can get this thing buffered?? Other software / workarounds / magic incantations / human sacrifice?? (And there'll be yet more steps after the buffer/dissolve, which would be no sweat with smaller data but will probably cause a minor earthquake if I try them with this beast).

EDIT :

After exploring this much more in-depth, I feel this question hasn't been adequately addressed by the duplicate question or the main suggestions here. Here's why:

The Dice tool appears to be for reducing the number of vertices in individual features? I don't have individual streamlines with >3,000 vertices, but in any case, I've used Dice with a vertices limit of 100, and still cannot run the dissolve.
Breaking the streams into 100 pieces and dissolving each one requires an "uber-Dissolve" at the end to combine the pieces, which will have buffers overlapping with one another. This "uber-Dissolve" experiences the same memory failure as the original approach, even though there is much less to dissolve after the 100 pieces have been processed. So this approach indeed doesn't seem to work, though for a completely different reason than what I stated earlier.

It seems PostGIS has implemented an approach called "cascaded unions" that addresses this general problem, which will be fun to explore, but in the meantime I was hoping for an ArcGIS-/Python-/QGIS-based solution.


texts: Sort by multiple fields and delete record with highest value

texts:
I have some GPS points that I am trying to find the highest precision point for. How would I go about sorting the table by time collected and horizontal precision and then deleting the highest precision for each record?

For example:

I would want to delete the record with Horz_Prec of 17.8 and keep the 13.8 one.


texts: QGIS problems getting layers to line up correctly.

texts: I am trying to add a raster layer and a vector layer to qgis.  I have set both of their CSR to NAD83 Oklahoma North State Plane.  It doesn't seem to work even with all the files set to those CSR.   What is likely the issue is that the raster has no extent, and QGIS doesn't seem to have a way to set the extent on rasters.  Coming from ArcGIS world where all these are set on the fly... Where else might I be going wrong.


texts: How to show Labels on Geoserver, only if they fit within the Polygon?

texts:
I'm trying to Set the Style via SLD on a Polygon Layer, such that it shows labels, and the Labels are visible only within the Polygon, and do not overflow outside the Polygon.

I've managed to get the Polygons to show, using the following SLD:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
  xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
    <NamedLayer>
        <Name>locality</Name>
        <UserStyle>
            <Title>locality</Title>
            <FeatureTypeStyle>
                <Rule>
                    <Title>locality</Title>
                    <PolygonSymbolizer>
                        <Stroke>
                            <CssParameter name="stroke">#47B8FF</CssParameter>
                            <CssParameter name="stroke-width">1</CssParameter>
                        </Stroke>
                    </PolygonSymbolizer>

                    <TextSymbolizer>
                        <Geometry>
                            <ogc:Function name="centroid">
                                <ogc:PropertyName>geom</ogc:PropertyName>
                            </ogc:Function>
                        </Geometry>
                        <Label>
                            <ogc:PropertyName>locality</ogc:PropertyName>
                        </Label>
                        <Font>
                            <CssParameter name="font-family">Arial</CssParameter>
                            <CssParameter name="font-size">11</CssParameter>
                            <CssParameter name="font-style">normal</CssParameter>
                            <CssParameter name="font-weight">bold</CssParameter>
                        </Font>
                        <LabelPlacement>
                            <PointPlacement>
                                <AnchorPoint>
                                    <AnchorPointX>0.5</AnchorPointX>
                                    <AnchorPointY>0.5</AnchorPointY>
                                </AnchorPoint>
                            </PointPlacement>
                        </LabelPlacement>
                        <Fill>
                            <CssParameter name="fill">#47B8FF</CssParameter>
                        </Fill>
                    </TextSymbolizer>
                </Rule>
            </FeatureTypeStyle>
        </UserStyle>
    </NamedLayer>
</StyledLayerDescriptor>

This has managed to show the Labels, but they still overlap adjacent polygons, like this:

How can I update this SLD, so that the Labels do not overlap adjacent polygons?


texts:
I have millions of seismic lines to publish on a map (3,000,000+) so WFS is not really an option at a large scale.  For this reason I will have to use WMS. My question is, is it possible using Geoserver as the engine and OpenLayers (latest versions) to present a map and allow the end user to click on an item and highlight the item?  Basic map functionality I would have thought.

There is a demo on the OpenLayers website but it appears this is the only example that does not work :( http://openlayers.org/dev/examples/getfeatureinfo-control.html

Thanks, Chris


texts: I am wanting to open shapefile files in a ASP.NET MVC application, was researching and found sharpmap but have not found examples of how to implement the MVC. Does anyone have any examples to help me develop.


texts: Marker appears in different position when zooming

texts:
I am using OpenLayers 4 and following the example http://openlayers.org/en/latest/examples/icon.html didn't work when I implemented. I need to place several maps from latitude and longitude coordinates.

For the moment I will leave as no zooming map, but I wish to know how to keep the mark in the same lat, lon position when zooming?

Position lost when zooming OpenLayers 4

My code:

        lat =  40.416775;
        lon = -3.703790;

        var iconFeature = new ol.Feature({
            geometry: new ol.geom.Point(ol.proj.fromLonLat([lon, lat]))
        });
        var iconStyle = new ol.style.Style({
            image: new ol.style.Icon(/** @type {olx.style.IconOptions} */ ({
              anchor: [0.5, 46],
              anchorXUnits: 'fraction',
              anchorYUnits: 'pixels',
              src: '/plugins/obs/gfx/marker.png'
            }))
        });
        iconFeature.setStyle(iconStyle);
        var vectorSource = new ol.source.Vector({
            features: [iconFeature]
        });
        var vectorLayer = new ol.layer.Vector({
            source: vectorSource
        });

        var map = new ol.Map({
            target: mapid,
            layers: [
              new ol.layer.Tile({
                source: new ol.source.OSM()
              }),
              vectorLayer
            ],
            view: new ol.View({
              center: ol.proj.fromLonLat([lon, lat]),
              zoom: 12,
            }),
            controls: ol.control.defaults({
               attributionOptions: {
                 collapsible: false
               }
            })
          });


texts: Is there any way to do (or emulate) the compositing operations that can be done in CartoCSS with CARTO VL? The closest I've been able to find (which hasn't been that close) is the ways to change the colour ramp and opacity, but I'm yet to work out a way to overlay/multiply/screen/etc the colours from two different layers in CARTO VL.


texts:
I'm taking a QGIS MOOC on Canvas, and in one of the labs we're to change layer style using an attribute table. However, the attribute table displays some very strange information on my system.  Attached is a picture... What's going on there?

Windows 7 Laptop - Acer
QGIS 2.6.1 Brighton

Attempts :

I have downloaded the data twice to replace it and started the lab from scratch.
I installed a package called "Table Manager" after reading SO, and restarted the software.
I changed the settings to "View all features."

At the very least, it appears I need the column headers to complete the lab assignment.

https://learn.canvas.net/courses/464/files/240563/download

You can download the data here.


texts: I want to mosaic a set of rasters using gdalwarp. I know I can first build a vrt and then use gdalwarp but, I have no idea if is it possible to do it without using vrt files? I also know about gdal_merge.py tool but I do not want to use this tool too. Is there any other method to do it?


texts: How do I catch PostGIS/PostgreSQL RAISE messages in QGIS?

texts:
I have a trigger set up in PostGIS/PostgreSQL which will show a RAISE WARNING when an inserted or updated polygon is within a certain distance of some predefined planning constraints. I want this warning to be sent to QGIS to alert a user to the proximity warning, but with a WARNING this does not occur. On the other hand, using RAISE EXCEPTION shows an error in the QGIS message bar but doesn't allow the insert/update to be committed (which I want - it should be a warning, not an error).

Example of the messages I want to be sent to QGIS:

CREATE TABLE input_table (
    id serial,
    geom geometry
             );

CREATE TABLE listed_buildings (
    id serial,
    geom geometry
             );

CREATE OR REPLACE FUNCTION check_constraint_proximity()
  RETURNS trigger AS
$BODY$
declare
dist_lb real;
begin
    dist_lb := (select st_distance(new.geom, st_collect(listed_buildings.geom)) from listed_buildings);

    if dist_lb <= 3 then
        RAISE WARNING 'Check planning boundary - it is within approx %m of a listed building', round(dist_lb);
    end if;
    if dist_lb <= 1 then
        RAISE EXCEPTION 'Check planning boundary - it is within approx %m of a listed building', round(dist_lb);
    end if;
    RETURN new;
end;
$BODY$
  LANGUAGE plpgsql;

CREATE TRIGGER check_constraint_proximity_trigger
  AFTER INSERT OR UPDATE
  ON input_table
  FOR EACH ROW
  EXECUTE PROCEDURE check_constraint_proximity();

insert into listed_buildings (geom) values (st_geomfromewkt('polygon((0 0, 0 1, 1 1, 1 0, 0 0))'));

insert into input_table (geom) values (st_geomfromewkt('polygon((2 2, 2 3, 3 3, 3 2, 2 2))'));
insert into input_table (geom) values (st_geomfromewkt('polygon((0 0, 0 1, 1 1, 1 0, 0 0))'));


texts: Load buffer by percentage QGIS plugin python error?

texts:
I couldn't install the buffer by percentage plugin on the QGIS 2.6 (I have Windows 8.1 on my pc), and gives a error message - the plugin is broken.

Then when I start the QGIS it gives the following window message:

Couldn't load plugin BufferByPercentage due an error when calling its classFactory() method

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\utils.py", line 208, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/LCP/.qgis2/python/plugins\BufferByPercentage\__init__.py", line 52, in classFactory
    from bufferbypercentage import BufferByPercentagePlugin
  File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\utils.py", line 460, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/LCP/.qgis2/python/plugins\BufferByPercentage\bufferbypercentage.py", line 38, in
    from processing.core.parameters import ParameterVector
  File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\utils.py", line 460, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named parameters

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.6.0-Brighton Brighton, exported

Python path: ['C:/PROGRA~1/QGISBR~1/apps/qgis/./python', u'C:/Users/LCP/.qgis2/python', u'C:/Users/LCP/.qgis2/python/plugins', 'C:/PROGRA~1/QGISBR~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\matplotlib-1.3.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\nose-1.3.3-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\tornado-4.0.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\backports.ssl_match_hostname-3.4.0.2-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\certifi-14.05.14-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISBR~1\\bin', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg']

I don't know anything about python.
Can anyone help me to fix this problem, so that I can use the buffer by percentage plugin?


texts: Create a colored hillshade in ArcGIS Desktop

texts:
Is it possible (and how) to create a colored hillshade using a color other than gray in ArcGIS Desktop as it is possible in ArcGIS Pro as described in this article by John Nelson on the ArcGIS Blog?
https://www.esri.com/arcgis-blog/products/arcgis-pro/mapping/kandy-kolored-multidirectional-streamlined-hillshade/


texts: Spatialite spatial queries do not function on GeoPackage

texts:
I have downloaded some data (as GeoPackages) from this excellent website: https://data.mfe.govt.nz/layer/52757-petroleum-wells/ and https://data.mfe.govt.nz/layer/52675-location-and-extent-of-nzs-aquifers-2015/ with the ultimate goal of writing a third GeoPackage where the two layers are related according to their spatial relationship (i.e. aquifers with references to intersecting wells).

QGIS is able to read these GeoPackages, and to display their geometries.

However when I use the Spatialite GUI (v4.3.0a) to query these datasets, I am running into odd geometry validation problems.

For example:

SELECT Count(*), GeometryType("geom"), Srid("geom"), CoordDimension("geom")
FROM "petroleum_wells"
GROUP BY 2, 3, 4

returns

1188    POINT    2193    NULL

i.e. that CoordDimension("geom") is NULL. This sensibly flows into literally every ST_* function I've tried: they all return -1 for any function operating on the geometries, which is explained in the documentation as: "–1 for UNKNOWN corresponding to a function invocation on NULL arguments." However, the geometry itself is not NULL. It is returned as the outcome of a simple select, albeit with a rather funny tpyoe: BLOB sz=47445 GeoPackegeGEOMETRY.

Oddly, ST_Is3D(geom) correctly returns 0 (false) when it could return -1 (unknown)... so it seems to have at least an understanding that's it's not in 3D.

ST_IsValidReason(geom) claims that the reason is Invalid: NULL Geometry even when I explicitly include a WHERE "geom" IS NOT NULL clause, and have successfully rendered the geometries in QGIS.

Is the GeoPackage geometry really invalid, or is there some issue with Spatialite, or some issue with the way I am using Spatialite?


texts:
I found two Raster Editors.

The first one is Serval Plugin and another one is  Editor Raster in SCP of QGIS. Editor Raster is working not well for me and Serval can edit only single pixel per time. However, I need to edit raster in a thousand of pixel.

So could you please recommend me some tool in QGIS that I can edit faster or kind of brush editor in QGIS?


texts: I am using ArcGIS 10.2 in order to create soil management zones for farmers. when I use the raster to polygon tool it creates small triangles out of the polygons and jagged edges that are very unrealistic looking zones. I have used the smooth polygon tool, as well as the nibble, focal statistics and majority filter while the files is still a raster in order to reduce the number of zones. Any ideas on an effective way to change these zones from raster to vector and look realistic?


texts: Combining three data extraction functions into one button?

texts:
I am looking to simplify my code. As of right now I have some code that allows users to extract GIS data in three different ways:

Extract by drawn extent
Extract by extent of the a county
Extract by extent by unit.

1st function:

function extractData() {
  //get clip layers
  var clipLayers = [];
  if (registry.byId("layer1").get("checked")) {
    clipLayers.push("SRA");
  }
  if (registry.byId("layer2").get("checked")) {
    clipLayers.push("Burn");
  }
  if (clipLayers.length === 0 || map.graphics.graphics.length === 0) {
    alert("Select layers to extract and draw an area of interest.");
    return;
  }

  var featureSet = new FeatureSet();
  var features = [];
  features.push(map.graphics.graphics[0]);
  featureSet.features = features;

  var params3 = {
    "Layers_to_Clip": clipLayers,
    "Area_of_Interest": featureSet,
    "Feature_Format": registry.byId("formatBox").get("value")
  };

  domStyle.set(loading2, "display", "inline-block");
  gp.submitJob(params3, completeCallback, statusCallback2, function (error) {
    alert(error);
    domStyle.set(loading2, "display", "none");
  });
}

2nd function:

function executeQueryTask(county) {
  var clipLayers = [];
  if (registry.byId("layer1").get("checked")) {
    clipLayers.push("SRA");
  }
  if (registry.byId("layer2").get("checked")) {
    clipLayers.push("Burn");
  }
  if (clipLayers.length === 0) {
    alert("No layers to clip, please select a layer you would like to clip...");
    return;
  }

  var county = document.getElementById("sel_county").value;
  var queryTask = new QueryTask("http://webgisdevint1/arcgis/rest/services/Alex_Try/Counties/MapServer/0");
  queryTask.on("complete", addToMap);
  var query = new Query();
  query.returnGeometry = true;
  query.outFields = ["NAME_PCASE"];
  query.where = "NAME_PCASE = '" + county + "'";
  query.outSpatialReference = map.spatialReference;
  queryTask.execute(query, function (featureSet) {
    var AOI = featureSet.features[0].geometry;
    var graphic = new Graphic(AOI, symbol);
    var features = [];
    features.push(graphic);
    var fSet = new FeatureSet();
    fSet.features = features;

    var params = {
      "Layers_to_Clip": clipLayers,
      "Area_of_Interest": fSet,
      "Feature_Format": registry.byId("formatBox").get("value")
    };

    domStyle.set(loading, "display", "inline-block");
    gp.submitJob(params, completeCallback, statusCallback, function (error) {
      alert(error);
      domStyle.set(loading, "display", "none");
    });
  });
}

And a third that looks exactly like the one above except that I query by units not county.

Each function is triggered by this:

registry.byId("extract1").on("click", executeQueryTask);

And this:

<button id="extract1" data-dojo-type="dijit/form/Button" data-dojo-props="iconClass:'Query', showLabel:false">
    Extract Data
</button>

Now here is my question: How can independently trigger each function from one and only one button?

As of right now I have three buttons that triggers each function.

Would it be possible to trigger each function independently from only one button?

For example, if users draw then the one and only "extract" button would extract data by the extent drawn by the user.

FYI I am using a dijit/from/dropdownfor selecting counties and units, which means there will always be one element selected in these two extraction processes.


texts: 1968 French Administrative Boundaries Shapefile

texts: I am looking for the 1968 French departments administrative boundaries. I've scoured the Internet but have had no luck.


texts: How to calculate the difference in area between two overlapping polygons?

texts:
I have two polygons that are overlapping each other. I need to calculate the area that's created in between these overlapping polygons. I used the "intersect" tool but it's not giving me the exact area. I used the field calculator function in the attribute table to manually subtract the areas between those two polygons. Is there a way to calculate the area between the overlapping polygons and represent it in a new polygon?

The area of the outer polygon is (red): 9.319671
The area of the inner polygon is (grey): 5.070906
The difference is: 4.248765 (need to represent this area in a new polygon)

The highlighted field in the attribute shows the value I manually calculated using the field calculator function.


texts:
This is related to an earlier question (Understanding U.S. Census MSA to place relationships?), on which I'm still unclear.
Can MSAs (CBSAs) split counties?

According to Glaeser et al. (Growth in Cities, 1992, p1135), in 1956:

  In some cases, an SMA* contains several counties; in others (only in
  New England), several SMAs split a single county.

*SMAs were MSAs of the time, and Glaeser et al. compare 1956 SMAs with 1987 MSAs.

But as of the latest U.S. census delineation file, I count (perhaps mistakenly) not more than one CBSA (Metro- or Micro-politan Area) per county, even in New England. This conforms with what the U.S. Census says, that "Counties or equivalent entities form the geographic 'building blocks' for metropolitan and micropolitan statistical areas throughout the United States and Puerto Rico." This could be because the U.S. Census's methodology has evolved on that score since 1956. For example, the latest methodology could prevent splitting counties by assigning ambiguous counties entirely to the 'dominant' CBSA (dominant in terms of commuting patterns).

Is it theoretically possible for contiguous CBSAs to split a county (especially an 'outlying county', per U.S. Census terminology)? This would also imply that contiguous CBSAs are possible.


texts: PostGIS: ST_Equals false when ST_Intersection=100% of geometry?

texts:
I have 2 datasets consisting of cadastral parcel data - roughly 125,000 rows each. The geometry column is WKB polygons representing parcel boundaries; all data are geometrically valid (the polygons are closed etc).

Some recent data arrived in a different projection to the base data being used for a comparison job - so I reprojected the newer one (base was 4326; the other was WGA94 that got brought into PostGIS as 900914... I reprojected it to 4326).

The first stage of the analysis was to find and store non-matching parcels; part of that is to identify and store parcels with identical geometry.

So I ran a very standard query (the code block below abstracts away schema details etc):

create table matchdata as
  select  a.*
  from gg2014 a, gg2013 b
  where ST_Equals(a.g1,b.g1)

ZERO results.

"Odd..." I thought. "Perhaps there have been tiny vertex shifts caused by the reprojection: that would be annoying, and really shouldn't happen."

Fortunately there is abundant aspatial data (5 identifier columns) that enable me to establish parcels that should be spatially identical: those with the same identifier, whose change-date in the 2014 table was before the max change-date in the 2013 data. That amounted to 120,086 distinct rows.

I stored the identifiers and geometries in a separate table (match_id), and ran the following query:

select apid,
       bpid,
       ST_Area(ag::geometry) as aa,
       ST_Area(bg::geometry) as ab,
       ST_Area(ST_Intersection(ag,bg)::geometry)/ST_Area(ag::geometry) as inta,
       ST_Area(ST_Intersection(ag,bg)::geometry)/ST_Area(ag::geometry) as intb
from match_id
order by inta

The first 16 values for inta and intb were identically zero, the next 456 were 0.99999999-ish (min 0.99999999999994, max 0.999999999999999), and rows 473 onwards were 1 - until row 120050, when the area of the intersection was greater than either geometry (the greatest value for inta and intb was 1.00000000000029, but still).

So here's my conundrum: if two geometries intersect spatially by between 99.999999999994% and 100.000000000029% of their respective areas, I would like "ST_Equals" to say "Yep.... I'll give you that one. Close enough".

After all, it's equivalent to being out by about 1 part in 16 trillion... i.e., as if the US national debt was off by less than 93 cents.

In the context of the circumference of the Earth (at ~40,000km), it's like being off by 0.0000000025km, tops (since to result in an area difference that small, any vertex shift must be even smaller).

According to TFD (which I have R'd) the tolerance for ST_Intersects() is notionally 0.00001m (1mm), so the implied changes in the vertices (which I confess I have not checked: I will ST_Dump() them and do so) would seem to be smaller than the tolerance. (I realise that ST_Intersects !== ST_Intersection(), but it's the only tolerance mentioned).

I have not been able to find out the corresponding tolerance for the vertex comparison undertaken by ST_Equals()... but it seems really odd that at least 120,000 of my rows ought to pass any sensible assessment of spatial identity, but don't.

(Note: I also did the same exercise using ::geography - with results that had more variability, but still more than 110,000 entries with a nice clean '1').

Is there a way to loosen the tolerance of ST_Equals, that doesn't require digging into the interstices of the code? I am not keen on doing that.

If not, is there a kludge that anyone is aware of?

Note: it would be good if the 'kludge' wasn't doing a bilateral comparison like

where ST_within(g1, ST_Buffer(g2, 0.0000001))
  and ST_within(g2, ST_Buffer(g1, 0.0000001))

   - I've done that: sure, it works... but it's a gigantic documentation PITA).

I can work around this, but writing the 20 pages to document the workaround - which will only ever come up again if we get dodgy data - is a PITA that I would rather not have to do given that it's likely to be a one-off.

(Versions: Postgresql 9.3.5; PostGIS 2.1.3)


texts: Buffer not calculated in CRS units despite reprojection

texts:
I have a point layer and a raster layer, where I'm trying to create a buffer around both of a certain radius in meters. I've tried saving each as the desired CRS (Albers Contiguous Equal Area) and used Warp for the raster layer. However, the buffer radius is still interpreted as degrees, as is the area of the buffer when I calculate that in field calculator. My current work-around is estimating degrees to meters, but this introduces error I'd prefer to eliminate. Is the problem the CRS I'm using, or am I not doing the reprojection properly?

Process for reprojecting:
Right click the layer in the layer tab > Export > Save features as > set CRS to Albers Contiguous Equal Area Conic, and save and add to project.


texts:
I'm building an application that stores locations with other data associated with those points. I'd like to be able to query an API which responds with a map with those points plotted out. (tangent) Eventually something like Leaflet will likely join the mix to explore those points dynamically.

Normally I'd use Google maps or Bing maps, but in this case we will be using custom maps so I'm leaning towards a hosting my own server to handle these requests. I'm leaning towards Geoserver or Mapserver. I'd prefer Geoserver but I don't seem to see any way to add points to a map I request, only the map itself. It seems like I could create points on a layer then add that, but then how do I create points/layers on the fly via the API...

Am I just missing something?

I'm open to suggestions if either of those two seem like a bad direction.


texts: XY_LINK field using Point_X & Point_Y Fields

texts:
In the below image, I am trying calculate my XY_LINK field to return the values from Point_X and Point_Y field.

Example of final result: {2474785.499662}{376075.063004}

My code is as follows:

Output(!POINT_X!, !POINT_Y!)

def Output(POINT_X, POINT_Y):
    FPX = round(float(X, 4)
    FPY = round(float(Y, 4)
    return "{%(FX)012.4f}{%(FY)012.4f}" % {'FX': FPX, 'FY': FPY}

The error I am receiving is as follows.  (I am brand new to Python)

Failed to execute. Parameters are not valid.
Error 000989: Python syntax error: Parsing error SyntaxError: invalid syntax (line 3)

I don't understand why the error is being returned on Line 3 and not like 2 also since they are the same with the exception of the X and the Y.


texts:
I have a raster that is tiled and I am looking to produce an average raster with a 3x3 window. I have tried three methods. The first method uses the standard ST_Mean4ma. An average raster is produced, however I can see the tiles in the result. This method took almost an hour to process. The second method I tried is a custom function that basically strips ST_Mean4ma to its base but obviously misses something as the tiles are even more visible - null values around the tile edges are produced. This method only took 23 mins. The third method I tried looked to ST_Union the raster first however it errored out with a memory allocation message.

How can I get a complete raster without tile edge effects and in a timely manner?

--TEST 1
CREATE TABLE p839.image_average AS
SELECT ST_MapAlgebra(ARRAY[ROW(rast, 1)]::rastbandarg[], 'st_mean4ma(double precision[][][], integer[][], text[])'::regprocedure, '32BF', null, null, 1,1
            ) as rast
    FROM p839.image_normal;
CREATE INDEX idx_image_average_rast ON p839.image_average USING GIST (ST_ConvexHull(rast));
--TIME = 51:40
--ARTIFACT TILES PRESENT = TRUE

--TEST 2
CREATE TABLE p839.image_average AS
SELECT ST_MapAlgebra(ARRAY[ROW(rast, 1)]::rastbandarg[], 'ImageAnalysis_Average(float[][][], text, text[])'::regprocedure, '32BF', null, null, 1,1
            ) as rast
    FROM p839.image_normal;
CREATE INDEX idx_image_average_rast ON p839.image_average USING GIST (ST_ConvexHull(rast));
--TIME = 23:53
--ARTIFACT TILES PRESENT = TRUE

--TEST 3
CREATE TABLE p839.image_average AS
SELECT ST_MapAlgebra(ARRAY[ROW(ST_Union(rast), 1)]::rastbandarg[], 'ImageAnalysis_Average(float[][][], text, text[])'::regprocedure, '32BF', null, null, 1,1
            ) as rast
    FROM p839.image_normal;
CREATE INDEX idx_image_average_rast ON p839.image_average USING GIST (ST_ConvexHull(rast));
--TIME = MEMORY ALLOCATION ERROR
--ARTIFACT TILES PRESENT = N/A

-- THE CUSTOM IMAGE AVERAGE FUNCTION
CREATE OR REPLACE FUNCTION imageanalysis_average(
    IN matrix double precision[],
    IN nodatamode text,
    VARIADIC args text[])
  RETURNS double precision AS
$BODY$
    DECLARE
    _matrix float[][];
        sum float;
    BEGIN
    _matrix := matrix;
        sum := 0;
        FOR x in array_lower(matrix, 1)..array_upper(matrix, 1) LOOP
            FOR y in array_lower(matrix, 2)..array_upper(matrix, 2) LOOP
                sum := sum + _matrix[1][x][y];
            END LOOP;
        END LOOP;
        RETURN (sum*1.0/(array_upper(matrix,1)*array_upper(matrix,2) ))::float ;
    END;
    $BODY$
  LANGUAGE plpgsql IMMUTABLE
  COST 100;

TEST 1:

TEST 2:

Anyone have anything on this?


texts: Unable to publish geocoding service [ERROR: 001369]

texts:
I manage to create a geo locator on arcmap and tested it.

but when I publish to arcgis server using arcmap or arcgis server manager. I encounter the error

1) Fail to publish service.

The error is very generic on logs.

may I know is there anyway to to troubeshoot this issue?


texts: QGIS Print layout: Making a Slider to create a number value

texts:
We currently have a QGIS template to deliver to a client,the thing is we are locking some items to make sure they are not edited. However some things like the Attachment number need to edited. Of course in theory you can adjust put an expression in a label, just like we are doing with the project number and project name. The thing is because the attachment number is always a number we were more interested in putting in a slider like item into the layout, like you have in the QGIS Forms if you define a range of numbers.
Is this possible or not? And if yes, how do I do this?
Edit: Ok so it is possible using this code more or less to get a slider in.
<!DOCTYPE html>

<meta name="viewport" content="width=device-width, initial-scale=1">

 Bijlage nummer:

    <input type="range" min="1" max="100" value="50" class="slider" id="myRange">
    <input type="range" min="1" max="100" value="50" class="slider" id="myRange">

The problem the slider is not interative in this HTML Frame, when I click on it I can only move the frame itself. How do I make it interactive?


texts: Converting a classification raster to a polygon, in order to calculate the area of each land use

texts: I have 6 land use categories. ultimately I would like to calculate the area of all of the classification categories in my layer e.g. the area of the deforestation, or the area of primary forest. I am trying to convert it from a raster to a polygon, in order to calculate the area from the attribute table. however when using the raster to polygon tool it is not working. I input the classified raster, leave VALUE as it is, then select an output in a new folder, tick SIMPLIFY POLYGON and 3 times it has not worked. I have looked at online tutorials and followed exactly what they do, and mine still does not work. Any suggestions please on what may be wrong?


texts: Visualizing layer using categorized AND graduated styles in QGIS?

texts:
I have a layer with road segments that come in five different classes and each have a value. I am trying to get a graduated color for the values (say, light blue for low and dark blue for high values) and at the same time make the highways (class 0) thicker than local roads (class 4) etc.

Now, usually I visualize the colors using about 10 classes (quantiles), so I'd end up defining somewhere around 50 rules: road class 1, value 0-100 - road class 1, value 100-500 and so on.

Is there a way to get there any quicker?

Using Version 2.18.10


texts: I've downloaded the GDAL-1.11.1 package but why can't I call on it in python?

texts:
I have a Mac and am very new to this coding stuff..

I'm using a code that includes

import osgeo as gdal

and so I've downloaded GDAL from https://pypi.python.org/pypi/GDAL/ which has a folder called "osgeo" and a python script called "gdal", but when I run the whole code I get the message "ImportError: No module named 'osgeo'"

Can someone give me the babysteps about how to import this correctly?


texts:
I have generated the map in google earth engine and QGIS. For example the 2019 NO2 levels for South Africa. I then filtered the layer for SA using a SA shapefile. I then try to export the layer. It then gives me Request payload size exceeds the limit: 4194304 bytes.

var band = NO2.select('NO2_column_number_density');
var banddate = band.filterDate('2019-04-01', '2019-04-30');
var median = banddate.median();
var SA = table.filter(ee.Filter.eq('locname', 'South Africa'));
var banddateSA = median.clipToCollection(SA);

var spectral = (['#a50026','#d73027','#f46d43','#fdae61','#fee08b','#ffffbf','#d9ef8b','#a6d96a','#66bd63','#1a9850','#006837']);

Map.setCenter(24.52, -28.77, 5);
Map.addLayer(banddateSA, {min: 0.0, max:0.0002, palette: spectral.reverse(Display)}, 'NO2');
Map.addLayer(SA, {color: 'yellow'});

Export.image.toDrive({
image: banddateSA,
description: "NO2",
fileNamePrefix: "NO2",
region: table,
scale: 30,
maxPixels: 1000000000000000,
fileDimensions: 10,
fileFormat: "GeoTIFF",
});


texts: What lens and camera settings are recomended for UAV (drone) aerial imagery capture?

texts: I have been toying with the idea of purchasing/bowrrowing a UAV and attempting to capture my own aerial imagery.  I have pretty good DSLR Canon gear and several lenses but I am wondering what lens, focal length and camera settings, including capture intervals, vehicle speed, etc. are recommended.


texts: ive created a point layer using the QChainage plugin and im able to change it it as i would a normal layer.. but the thing is that after i save and exit and re enter the qgis file the layer is gone .. is there a way to save it as an shp file so i can load it later on and use it with different files?


texts:
I have a fairly simple question. What is the best strategy to do the following scenario? I think this is a very common scenario:

I have a several thousands of dataset with polygon/points. Each dataset/record can have multiple polygons.
Allow users to draw a polygon and search which datasets fall into that space
I have several predefined polygon/areas, that the users can also use/pick to search for datasets.

When the search is done, just show the list of the datasets, display as markers on a map. I'm not worried about this bit because it has been done just with OpenLayers.

My current options (I am new and some of my current options may not even be possible. Please let me know if that is the case) are:

Put datasets in PostGIS, set it up in Geoserver. Put predefined polygons straight in geoserver. Use WMS filtering.
Put datasets in PostGIS, but not in Geoserver. Put predefined polygons in PostGIS and Geoserver for displaying. Use PostGIS intersection.
I've heard from another person that for performance ,it is best to calculate the intersection using a custom app, e.g. Java / other

I'm not sure what is the best in terms of performance and efficiency. Suggestions are very welcomed!


texts: Get PostGIS geometry viewer route drawing

texts:
I am using pgrouting pgr_dijkstra function to calculate the shortest path between two points. Sorting the result by "seq" gives me the correct sequence of the path to my destination. In the result I can also see the graphical output of my query directly on pgAdmin, with the help of the st_length (geom_way) function. With the geometry viewer, a route is drawn, as if it were a Google Maps.

Through the development of a program in Java or any other language, is there any way to download this image with the mapped route or is there an API that I can pass the data returned by the pgr_dijkstra function and that returns an image with the path to insert into a program?
I just needed a static image, a "print" of that route.


texts: Assigning multiple Customize keyboard shortcut keys for ArcMap

texts:
I am using ArcMap 10.1 and 10.2.1.
I have attached the below image for reference.  I want to assign multiple Customize keyboard shortcut keys at the same time.
Is it possible do this using ArcPy/Python?
This is to save effort changing the Customize keyboard shortcut keys from my MXD file for other people.


texts: ArcGIS and Excel - How to get Text format to stick when importing?

texts:
I have a table with data based on intersections of a cities general plan, zoning, and parcel data. I exported the resulting information into a table and cleaned it up in Excel. Now I want to bring it back into ArcMap and table join this data back into a clean parcels feature class using the APN numbers. The problem I am having is that the APN column includes a string of 8 numbers AND letters which, when brought into Arcmap, leaves the number-letter APN's as NULL.

Likely, this is because bringing this table back into Arcmap makes the column LONG when it should be TEXT to maintain the format of number-letter for purposes of table joining.

I can not, for the life of me, get the 'text' cell format in excel to stick when I bring it back into ArcMap. It always changes it to LONG which then results in NULL values.

I have even tried to save it as .txt rather than .csv and this still happened.


texts: Suggestion for java based GIS sdk for desktop application

texts:
I am developing a desktop GIS application using JavaFX that needs to run on Linux. I have been using ArcGIS java sdk and it is really good for all my requirements. But when I started testing it on our actual field machines it is not working. My client is using an AMD driver on their Centos machines and they not ready to change it because they say they found this as a solution to many of their problems. Some how the ArcGIS java sdk not working very well with the AMD driver so I am searching for an alternative. I really liked ArcGIS it has all the features that I required. So far I was not able to find a good alternative for that, can somebody pleas suggest an alternative solution.

For reference on of my coworker posted the problem in another thread https://gis.stackexchange.com/questions/114488/problem-with-arcgisruntimesdk-for-java-on-centos-6-2-after-installing-amd-cataly

I tried with the latest 10.2.4 version still not working.


texts:
In QGIS I have a rule-based label. I ALSO want to set the blocking (prevent the label from being overlapped by other features / labels). However, when I try to do this, my rule-based labeling is disappeared.
Is this a bug or a feature?


texts: How can you tell which port number an existing Database Connection is configured for?

texts:
I have a database connection to a postgresql database on a non-default port, 5433.

I specify TESTDB,5433 as the instance and I can connect with no errors.

If I open that connection properties again, the Instance is displayed with no port configured:

When I hit cancel and the double-click the connection, I connect successfully to the 5433 port (confirmed by checking the server status of the 5433 and seeing the incoming connection). Even after hitting OK on the Instance without a port...the connection still connects to 5433.  Only after expliciting specifying TESTDB,5432 can I switch back to the 5432.

Since both instances have the same data schema (and often the same data), how would an ArcGIS user know which instance he or she is connected to?


texts:
I am getting exactly this error when trying to add locations to a service area within ModelBuilder.  However, the symptoms mentioned in that link relate to my tool being called from a differing language, and mines the standard US install.

I have the iterator pulling a table in (with a date field) and pushing it to a 'time of day' parameter for the Make Service Area.  I think thats all setup fine.  Just not sure why its having issues with adding facilities (facility is one point from a featureclass) and when adding the previous layer to the display, I can see there is a sublayer called Facilities.

If I remove the Iterator, it works.
With the iterate field values, it works if I don't use an inline variable for changing the output layer name in make service area to Service Area%i% (i.e it overwrites the output every time).

Not too sure where to go from here.

Can anyone provide me with some tips on how to get this working?


texts: How to configure a custom ArcGIS Pro Add-in

texts:
How can ArcGIS Pro add-ins be configured?

I am building an add-in for ArcGIS Pro. This add-in will access external resources via HTTP calls to a SOAP service i have built. I want to be able to "point" the add-in to a specific version of the SOAP service for testing. And then i want to be able to change that configuration.

How can i do this with my add-in?

I would like to avoid having to hardcode the url to the external resource into my built add-in file.


texts: Leaflet - create accordion based on GeoJSON data

texts:
My webapplication displays a map and a sidebar. Within the sidebar I would like to have an accordion that displays the 10 items that are also displayed as points in GeoJSON format on the map.

For the html page I have

and the function initDescription() is called when the body is onload.

The function initDescription() accesses the JSON file and builds the accordion elements.

for (var i=1; i<=Object.keys(descriptionJsonFile).length; i++){

document.getElementById("description").innerHTML +=
    "<button class='accordion'>Trail Number" + wegbeschreibung[i].Titel; + "</button>" + "
" + "
" + description[i].TrailDescription+ "

" + "
";

Finally, in the main JavaScript file it is handled what happens when the accordion element is clicked.

var acc = document.getElementsByClassName('accordion');
var panel = document.getElementsByClassName('accordion-panel');

for (var i = 0; i<acc.length; i++) {
acc[i].onclick = function() {
    console.log("accordiuon is clicked");
    var setClasses = !this.classList.contains('active');
    setClass(acc, 'active', 'remove');
    setClass(panel, 'show', 'remove');

    if (setClasses) {
        this.classList.toggle('active');
        this.nextElementSibling.classList.toggle('show');
       }
     }
   }
   function setClass(els, className, fnName) {
    for (var i = 0; i < els.length; i++) {
        els[i].classList[fnName](className);
    }
}

The data from the JSON file are read out correctly and the accordion elements displayed, but it does not work to click and expand the accordion elements. I already tried to get the information via document.getElementById("description").getElementsByClassName('accordion') but with no success.

Anyone how has an idea what goes wrong or what I could try?

EDIT:
I found a workaround: putting the function initAccordion() directly as script within the html page at the desired position and not as an external script file, it works.

Anyone who knows why the body onload is called after the other content in the html body is loaded? Does the script directly within the html need more runtime than calling it externally?


texts: How to automatically transform between base layer and vector layer in OpenLayers

texts:
I have a simple example on http://jsfiddle.net/29bsw/

I have a map with an OSM layer as base layer and an additional vector layer. The projection of the map is defined as EPSG:900913 and the projection of the vector layer is defined as EPSG:4326:

//Map
var map = new OpenLayers.Map("map", {
    projection: new OpenLayers.Projection("EPSG:900913")
});

//Base Layer
var osm = new OpenLayers.Layer.OSM();
map.addLayer(osm);

//Vector Layer
var vectorLayer = new OpenLayers.Layer.Vector("Simple Geometry", {
    projection: new OpenLayers.Projection("EPSG:4326")
});
map.addLayer(vectorLayer);

Now I would like to add some features to the vector layer without transforming any geometry point from EPSG:4326 to EPSG:900913. The point should be lie in Hamburg, Germany, but without transforming (see commented on line 21 in the example) the point lies near by Africa.

Because I specified the projection of the map and of the vector layer, is there any way that OpenLayers transform the point from EPSG:4326 to EPSG:900913 automatically?


texts: Removing selection area from canvas when no layer is available or selected

texts: I am new to QGIS 3 and I can not figure out the way to remove the selection from empty canvas or with Google Maps layers selected this area can not be removed. I tried to remove it by selecting a layer but still not successful.


texts: Coordinates do not export in UTM format I can use

texts: I am trying to export my attribute table so I can view the UTM coordinates, which I need to present in a table in my report. I have got as far as exporting the X and Y coordinates, using the same projection as my project, but the results come in as long strings of numbers that I cannot find information on how to convert. I would really like the output to be as zone, easting, northing (Zone 10U). I have attached the export window and the CSV table result.


texts: Why is shapefile not recognised by Mason?

texts:
I have a shapefile of an Ordnance Survey Mastermap ITN road network that I am trying to load into an agent-based modelling software called Mason. When I try to run my model I get the following error:

  Exception in thread "main" java.lang.IllegalArgumentException: shpFile
  is null; likely file not found

My model works fine with every other road network (OS Mastermap Meridian 2, Boundary-Line etc.) I can throw at it and can display other shapefiles I'm using. All the shapefiles show up in QGIS just fine.

The OS ITN data was downloaded as a compressed .gz file, with two GML files within. These were were then converted to .MIF/.MID using OSM2MIF. These were then loaded into QGIS and saved as .SHP (right click>Save As...>ESRI Shapefile).

In QGIS and ArcMap the resultant shapefile of the road network looks as expected and all the metadata is present. However, when I try to load it using my Mason model, it fails. It seems to only fail using this shapefile so I know it's not my model and is thus a problem with the shapefile.

Have I done something wrong during the conversion? Does anyone know how best to use Ordnance Survey Mastermap ITN and convert to usable shapefiles?


texts: Generate point of view (POV) image inside Digital Surface Model (DSM) or 3D City model

texts:
I'm currently working with DSMs (Digital Surface Models) and user perspective inside that model.

I'm looking for a Python library, which lets me generate a specific view (from the POV perspective) depending on coordinates and orientation.

Let's assume the following conditions:

The user and model have the same coordination frame and I'd like to know what the user sees when he's standing on this model, but I don't want to render a 3D model of that perspective, but a 2D view.

It's like the user would take a photo from that point.

The output product should look like this:


texts:
I have two shapefiles with the same projection, but they could not overlay each other. I understand that the extent is different, but I could not figure out how to change it.

Here are the two layers properties:


texts: Batch attribute editing across multiple layers

texts:
Using ArcMap 10.5

I want to find a way to do a bulk attribute update across multiple layers.

Example:

I have 10 layers in a Storm Network and each layer has 5 identical fields (project related information) that I would like to edit all at once. Meaning I enter the data for each field once, and a tool or simple script applies the data to every selected asset across all layers.

My hope was that there already was an add-in somewhere out there or a script that could be downloaded and configured. But I haven't been able to find anything like this which I find surprising.

Anyone willing to share their process for achieving bulk updating of attributes?


texts: Keeping symbology style in a polygon shapefile when converting to KML

texts: I have a polygon shape with 23 classes and have created patterned symbology styles which I will like to have retained when I open in Google earth pro. Is there any how I can achieve this in either ArcMap or QGIS. My programming skills are very poor.


texts:
I am new to OpenLayers. I need guidance on below.

I am able to create map using OpenLayers and able to create polygon interactions on the map and able to extract as GeoJSON object with coordinates.

Now I need to send this GeoJSON request to my GeoServer to get the intersected grids on the OpenLayers (grids are configured in GeoServer).

How can I achieve this using RESTful API?


texts: Extracting specific values to new raster in ArcGIS Desktop?

texts:
I have a map with data stored in file geodatabase (gdb) format. As I opened it in ArcGIS 10, I see that it has data for 12 years, each year encoded in different color.

How do I create 12 different maps from this one big dataset, when in the attribute table there is only an objectid column, a value column (12 levels, corresponds to each year), and the count number column (the number of cells that has data for each year)?


texts: Linestring table (of MultiLinestrings) ST_Intersection

texts:
I'm moving from GeoPandas to PostGIS and GeoAlchemy.
I've a table full of multilinestring geometries (timezones from NaturalEarth.com) and I've a linestring which represents what will be the route of an aircraft. I would like to get the intersected timezones AND the position of intersections.
So far I have the intersected timezones using
class TimeZoneLine(Base):
    __tablename__ = 'TimeZone (LineGridExp)'
    __table_args__ = {'autoload': True}

lineGeom = func.ST_GeomFromText('LINESTRING(-14 44, 6 44)',4326)

query = session.query(TimeZoneLine).filter(func.ST_Intersects(TimeZoneLine.geom, lineGeom))

for row in query:
    print (row.zone)

which prints
<__main__.TimeZoneLine object at 0x107656dd8>
<__main__.TimeZoneLine object at 0x107656e48>
<__main__.TimeZoneLine object at 0x107656f98>
<__main__.TimeZoneLine object at 0x107683128>
<__main__.TimeZoneLine object at 0x107683278>
<__main__.TimeZoneLine object at 0x1076833c8>

which are the timezones intersected. My question is, how to get the Intersection points.
This seems to have some thing similar but the below is closest I've come:
query =  session.query(select([column('zone'),column('id'), column('iso_8601'),column( 'map_color6'), func.ST_intersection(lineGeom, TimeZoneLine.geom)]).select_from(TimeZoneLine).where(func.ST_Intersects(lineGeom, TimeZoneLine.geom)))

but this returns
sqlalchemy.exc.ProgrammingError: (psycopg2.errors.SyntaxError) subquery in FROM must have an alias


texts: Why does Blender not accept DT2 file (elevation data from EarthExplorer)?

texts:
Trying to do the Creating Shaded Relief in Blender tutorial

https://somethingaboutmaps.wordpress.com/2017/11/16/creating-shaded-relief-in-blender/

Got the file from SRTM 1 Arc-Second Global. Got also a DEM file from NASA STRM3 STRMGL30 and Blender does not read it. The file does not appear in Blender's file opener.


texts: Correct way to delete a field, add a field and update attributes in pyqgis

texts:
I am working on a script that runs a model and then appends the result to a spatialite layer. The model runs zonal statistics using 4 different rasters and appends them to the attribute table of a vector polygon layer with a single feature. The attribute table of the model output (habitatModel) looks as follows:

id  | hs90_mean | hs14_mean | ls90_mean | ls14_mean
----+-----------+-----------+-----------+----------
 1  | 46.043816 | 44.009554 | 47.859847 | 45.308090

The type of the id field is integer, and I need to replace that with a string identifier, so I load the layer as a QgsVectorLayer, add a genspec field and delete the id field:

habitatModelLayer = QgsVectorLayer(habitatModel,genspec,'ogr')
habitatModelLayer.dataProvider().addAttributes([QgsField("genspec", QVariant.String)])
habitatModelLayer.dataProvider().deleteAttributes([habitatModelLayer.fields().lookupField('id')])
habitatModelLayer.updateFields()

Then I try to update the attribute values using dataprovider().changeAttributeValues():

iterator = habitatModelLayer.getFeatures(QgsFeatureRequest().setFilterFid(0))
feature = next(iterator)
attrs = {habitatModelLayer.fields().lookupField('genspec'):'2759-9',
    habitatModelLayer.fields().lookupField('hs90_mean'):feature['hs90_mean'],
    habitatModelLayer.fields().lookupField('ls90_mean'):feature['ls90_mean'],
    habitatModelLayer.fields().lookupField('hs14_mean'):feature['hs14_mean'],
    habitatModelLayer.fields().lookupField('ls14_mean'):feature['ls14_mean']}
habitatModelLayer.dataProvider().changeAttributeValues({0:attrs})
habitatModelLayer.commitChanges()

I then add the layer to the spatialite file database using qgis:importintospatialite:

processing.run("qgis:importintospatialite", {
    'INPUT':habitatModelLayer,
    'DATABASE':spatialiteLayer,
    'TABLENAME':'models',
    'PRIMARY_KEY':'genspec',
    'GEOMETRY_COLUMN':'geom',
    'ENCODING':'UTF-8',
    'OVERWRITE':True,
    'CREATEINDEX':True,
    'LOWERCASE_NAMES':False,
    'DROP_STRING_LENGTH':False,
    'FORCE_SINGLEPART':False
})

The resulting layer has the correct values for the zonal stats, but the genspec field is empty. Any idea why this would be?

A strange side-note is that when the habitatModel layer is added to the canvas as a result of the script, it initially still shows the new attributes squashed into the old fields (and trying to conform to the old types):

id  | hs90_mean | hs14_mean | ls90_mean | ls14_mean
----+-----------+-----------+-----------+----------
 46 | 44.009554 | 47.859847 | 45.308090 | 2795

but after a few minutes the attribute table corrects itself:

 hs90_mean | hs14_mean | ls90_mean | ls14_mean | genspec
-----------+-----------+-----------+-----------+---------
 46.043816 | 44.009554 | 47.859847 | 45.308090 | 2759-9

Is there some kind of buffer that needs to be flushed before I write to the spatialite database?


texts:
I am trying to construct a geodetic line between destinations. I have followed two guides without luck. This are the steps I've done:

I created a polyline shapefile and named it "Distance"
I went to advanced editing a pressed "construct geodetic"
Now there should be a template called "Distance" to choose but as you can see in the picture it isn't there.

What am i doing wrong?


texts:
The Mapillary API for sequences https://www.mapillary.com/developer/api-documentation/#the-sequence-object doesn't make any mention of street names. This seems strange for a street images service.
Nevertheless, I'm trying to get a set of Mapillary sequences for any given route, say I-75 South. I could go to the viewer and find sequences and try to patch together some from a given user. Or,  I could use the API to look for sequences in a bounding box and then try to stitch those together. Is there any better way?


texts:
I am trying to create a graph with the cloud cover over time, per year. The following code gives me a graph, but I don't manage to get the dates on the x-axis line up chronological. Am I doing something wrong?

Secondly, this line contains 2 years of Landsat 8 data, I would like compare the annual values with the averaged values of the entire time series (in this case 2 years). Any ideas how to do so?

// ------CLOUD COVER CHART--------------------------------------------

var l8toa = l8.filterBounds(aoi);
var results = l8toa.map(getCloudScores);
var chart = ui.Chart.feature.byFeature(results)
    .setChartType('LineChart')
    .setOptions({
      title: 'Cloud Cover over Time',
      hAxis: {
        title: 'Date',
 /*      ticks: [{v: 1, f: 'January'},
                {v: 4, f: 'April'},
                {v: 7, f: 'July'},
                {v: 10, f: 'October'}]*/},
      vAxis: {
        title: 'Cloud Cover'
      },
      lineWidth: 1,
      pointSize: 3
});


texts: Why is the CORS coverage density in the US is so high?

texts:
According to the data I found online, it appears that the density of CORS (Continuously Operating Reference Stations) in the US is significantly higher compared to other countries. Why?

(NOAA coverage map)

(EUREF coverage map)

(Australian regional map)

(IGS network)


texts: Python parser in field calculator malfunctioning

texts:
Shouldn't this return everything but the last word in the street name?

!Street_Name!.split(" ")[:-1]

Python parser is checked. Receiving nulls after calculation. If I try to slice[0], I get the first word. If I try to slice [-1] I get the last word. Why won't [:-1] return anything?


texts: Finding nearest points based on distance and the user location and display on a Folium or Leaflet map?

texts:
I want to display the nearest 'services' points using the distance function on the Folium map. I need to pass the user location from the form and find the nearest services based on the current location. This is what I have my views.py file:
from django.shortcuts import render, get_object_or_404
from .forms import addressModelForm
from django.views import generic
from django.views.generic import TemplateView
from django.views.generic.list import ListView
from django.core.serializers import serialize
from django.http import HttpResponse
from django.contrib.gis.geos import Point
from django.contrib.gis.db.models.functions import Distance
from .models import Services, GautengPlaces, Servicesgp, hospitals, Pointdistance
from geopy.geocoders import Nominatim
from django.contrib.gis import geos, measure
from django.contrib.gis.db import models
import folium
from urllib import request

# Parktown long/lat
latitude = -26.167
longitude = 28.033

# Turffontein long/lat
#latitude = -26.2446
#longitude = 28.0397

user_location_init = Point(latitude, longitude, srid=4326)

# Create your views here.

def calculate_distance_view(request):
    form = addressModelForm(request.POST or None)
    glocator = Nominatim(user_agent='reporter')
    model = Services
    context_object_name = "services"
    m = folium.Map(width=800, height=500, location=user_location_init, zoom_start=10)

    if form.is_valid():
         user_address = form.cleaned_data.get('user_location')
         user_location = glocator.geocode(user_address)
         d_lat = user_location.latitude
         d_long = user_location.longitude
         current_location = Point(d_lat, d_long, srid=4326)
         nearbyobj = Services.objects.annotate(distance=Distance("location", current_location)).order_by("distance")[0:3]

         m = folium.Map(width=800, height=500, location=user_location_init, zoom_start=10)
         folium.Marker([d_lat, d_long], tooltip='click here for more', popup=user_location,
                       icon=folium.Icon(color='red', icon='cloud')).add_to(m)

    m = m._repr_html_()

    context = {
        'form': form,
        'model': model,
        'nearbyobj': 'nearbyobj',
        'map': m
    }
    return render(request, 'index.html', context)

Models.py to pull the services from the model
class Services(models.Model):
    name = models.CharField(max_length=254)
    location = models.PointField(srid=4326)

    def __unicode__(self):
            return self.name

    class Meta:
        verbose_name_plural = " Services"

I'm not sure how to go about using 'nearbyobj = Services.objects.annotate(distance=Distance("location", current_location)).order_by("distance")[0:3] '  to display the nearest services based on the value we get from the form(user_location). I'm fairly new to GeoDjango, trying to get this project up and running.


texts: Convert KML to shapefile by ArcGIS ModelBuilder?

texts:
I am trying my hand at ModelBuilder today.  This is a simple task; I wish to accomplish this to:

Make this task more autonomous
Hone my ModelBuilder Skills.

All I'm trying to do is create a simple tool that will allow me to convert a KML/KMZ (or several via batch) to a shapefile which will subsequently allow me  to save the shapefile in the desired folder.  When I attempt to do this process manually I have to convert via 'KML to layer' tool then export layer to shapefile by right clicking "export data" in the TOC and then adding it to the layer.  This process becomes tedious if I have several KMZs.

Can anyone assist?

I am using ArcGIS Desktop 10.5.

Below is my first attempt which failed miserably.


texts: Why .getOutput(0) is needed after GetCount_management?

texts:
parcelCount = int(arcpy.GetCount_management("Parcels_lyr").getOutput(0))
print parcelCount
while parcelCount < 200:
    arcpy.SelectLayerByLocation_management("Parcels_lyr","WITHIN_A_DISTANCE",\
                                           "Parcels_lyr","150","ADD_TO_SELECTION")
    parcelCount = int(arcpy.GetCount_management("Parcels_lyr").getOutput(0))
    print parcelCount

The above was part of code. I don't understand the use and the need of getOutput(0). A reference said "because parcelCount is a python object and not a simple variable, you must add .getOutput(0) at the end of the object to reference the value's index number." I am totally confused. I look up online, including python doc, but seems not so much helpful info available.

I also searched this site. There is 1 post mentioning it, but not much.

Can anybody explain it?


texts:
I'm trying to find a good source of NDVI imagery, but I'm struggling to find proper high res imagery, particularly for northern Europe.

Does anyone have references to high res imagery, that is frequently updated?

I've tried looking at Landsat 8 imagery, using the EarthExporer, but there's no imagery in the Vegetation datasets for Europe.

I'd like to find imagery in the 50-100m/px resolution, but I don't know if that's possible without being charged an arm and a leg for it?

If no pre-processed NDVI imagery is available, it would be fine to get raw imagery that I could use to process NDVI imagery from.


texts: Reclassify raster pixels underlying certain polygons using Python

texts:
I've got a shapefile with many polygons, each of which exactly overlays a contiguous region of raster cells; an example in yellow:

For some of these, I want to use the polygon as a sort of analysis mask, so that I can reclassify only those cells that are overlaid by the polygon according to one of the polygon's attributes. I have to do this thousands of times to the same raster, so I'm looking for a Python-automated way to do it, using (hopefully) any of GDAL/OGR, SAGA, QGIS, or numpy arrays.


texts: Graticule will not display in layout view

texts:
Problem:

I am trying to create a graticule for a map in ArcMap 10.7.1. When I create the graticule and hit "apply", the graticule does not display in layout view or when I export the map.

Steps to create the graticule:

View -> Data frame properties -> Grids
New Grid -> Graticule -> Labels Only -> 5 degree meridians / parallels
Etc.

Attempts to resolve:

Changing the Label Offset point value to at least half the Label Size (from here).
Defining a coordinate system for the data frame (from here).
Turn on grid (View -> Grid)

I have tried reducing intervals to minutes or seconds, but graticules are still not appearing.


texts: Creating network datasets from OSM data in QGIS?

texts:
I am having a difficulty to use network analysis in ArcGIS because of unavailable network datasets.
So, I googled it and realized that it is possible by using OpenStreetMap.
However I don't know how to create a network dataset from OSM in qgis and then incorporate it to ArcGIS.
I am a very beginner so please let me know in details.


texts: Where water will be on a mudflat as tide comes in/goes out?

texts: My study area is large intertidal mudflat in NE England. I want know where the water will be depending on tide height. The attached shows an aerial of my study area. I've photoshopped in the blue lines to show the type of information I need, therefore where the water will be depending on water height as the tide comes in/goes out. Is this data available anywhere?