texts: Methods to more efficiently generate regularly spaced research points inside of a multipolygon

texts:
I'm currently using the tried-and-true, but extremely slow method of:

via QGis, generating two (ESRI shapefiles) of a large number of research points over a multipolygon which are generated using a bounding box of the extremities of the multipolygon (a huge square shaped matrix of millions of points). This doesn't take that long.
Then (also in QGis) running an intersection on the "raw" layer of research points on the actual multipolygon areas to output a layer of research points that only lie within multipolygon areas. Takes a very, very long time and seems very inefficient.

This still works well enough using sufficiently spaced out research points, but I'm now running an intersection of tens of millions of research points to output intersecting points totaling perhaps only 10,000. Its going to take over 8 hours* on my 16 core CPU. Is there a better way to do this (an optimized intersection operation) that won't cause too much brain damage?

*edit: took nearly 20 hours, and I know there's no way I used more than one core of my fancy processor. Arghh


texts: Forcing ArcMap to display HTML pop-ups for just one layer?

texts:
At 10.5.1, I have an .mxd that contains several polygon feature class layers.  I have created an HTML pop-up (my first-ever!) for one of those layers.  Let's call that layer "A".

The pop-up works great, in that it points to the correct url.

But there's a problem: when I click on a feature, the resulting pop-up displays data concerning the other layers - not the desired layer "A". I can force layer "A" to pop-up only if I turn off the display for all the other layers, leaving only layer "A" checked on in the TOC.

Turning the other layers on and off is awkward and time-consuming.  How do I specify that the pop-up tool only apply to layer "A" in the TOC?


texts: Getting points digitized in Web AppBuilder geoprocessing widget back over to original coordinate system for calculations?

texts:
This is a follow-up question to Having two input feature layers and derived output feature layer in ArcGIS Python Toolbox tool?.

The following Python Toolbox tool accepts one or more points from the user and calculates the station, offset, and orientation of the point(s) based on a single polyline (an 'alignment').  It works fine as a geoprocessing tool in ArcMap with all the layers in the local Tennessee State Plane NAD83 coordinate system.  However, I'd like this to also work correctly from within the Geoprocessing Widget in Web AppBuilder.  When the user goes to digitize the point(s) in WAB, they are of course in the Web Mercator (auxiliary sphere) coordinate system.  For the geoprocessing service to work correctly, I need to somehow project those points from Web Mercator back to TN State Plane so that the queryPointAndDistance method correctly calculates the station, offset, and orientation. What is the best (or easiest) method for accomplishing this?

import arcpy

class Toolbox(object):
    def __init__(self):
        self.label = "Station Offset toolbox"
        self.alias = "StationOffset"

        # List of tool classes associated with this toolbox
        self.tools = [CalculateStationOffset]

class CalculateStationOffset(object):
    def __init__(self):
        self.label = "Calculate Station Offset"
        self.description = "Calculate Station Offset"

def getParameterInfo(self):
    #Define parameter definitions

    # Input Features parameter
    in_features = arcpy.Parameter(
        displayName="Input Features",
        name="in_features",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")

    in_features.filter.list = ["Point"]

    # Alignment parameter
    in_alignment = arcpy.Parameter(
        displayName="Alignment",
        name="in_alignment",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")

    in_alignment.filter.list = ["Polyline"]

    # Derived Output Features parameter
    out_features = arcpy.Parameter(
        displayName="Output Features",
        name="out_features",
        datatype="GPFeatureLayer",
        parameterType="Derived",
        direction="Output")

    out_features.parameterDependencies = [in_features.name]
    out_features.schema.clone = True

    parameters = [in_features, in_alignment, out_features]

    return parameters

    def isLicensed(self):
        return True

    def updateParameters(self, parameters):
        if parameters[0].altered:
            parameters[1].value = arcpy.ValidateFieldName(parameters[1].value, parameters[0].value)
        return

    def updateMessages(self, parameters):
        return

def execute(self, parameters, messages):
    fc_pnt         = parameters[0].valueAsText
    fc_line        = parameters[1].valueAsText

    # Fields
    fields = ['SHAPE@', 'MEAS', 'Distance', 'NEAR_X', 'NEAR_Y', 'Orient', 'POINT_X', 'POINT_Y', 'NEAR_ANGLE']

    # Get line geometry - assumes only one feature in feature class
    polyline = arcpy.da.SearchCursor(fc_line, "SHAPE@").next()[0]

    # Loop over the point feature class
    with arcpy.da.UpdateCursor(fc_pnt, fields) as cursor:
        for row in cursor:
            dist=polyline.queryPointAndDistance(row[0], False)
            row[1] = dist[1]
            row[2] = dist[2]
            row[3] = dist[0].centroid.X
            row[4] = dist[0].centroid.Y
            if dist[3] == 0:
                row[5]="Left"
            else:
                row[5]="Right"
            row[6] = row[0].centroid.X
            row[7] = row[0].centroid.Y
            print(dist)
            cursor.updateRow(row)


texts: Using raw_input python to accept string with dash (-)?

texts:
I have the following python input. The user will enter a value with a dash always. E.g. 2015-0001

Var1 = raw_input("Enter the number with dash: ")
arcpy.CalculateField_management("all_reports.gdb\\joined", "YearNo", Var1, "PYTHON_9.3", "")

"YearNo" is a text field. If I do "print type(Var1)" it tells me that it is a string.

When I look at my YearNo field is shows 1972. It is doing a subtraction of the numbers rather than keeping a string together.

Is there a way around this?


texts: Openlayers 6, App crashes when I set too many text labels

texts:
I'm looking for an efficient way to set conditionally text for over 3000 features. Below is a little code on how I do it now, but my Android App freezes and crashes sometimes because it exceeds the memory limit of 2 GB. I did try to cache my styles first by color but I end up with only two styles for every feature so I can't set text individually. I did notice that, when I set the text font size to 24px, my app crashes on start, when I return the style function, because of that I set the font size to 18 px.

var layer = new ol.layer.VectorImage({
            source: vectorSource,
            style: function(feature) {
                style.getFill().setColor(getColorByID(feature.get('id'), feature.get('pe_nr')));
                return style;
            }
});

var style = new ol.style.Style({
  fill: new ol.style.Fill({
    color: '#FFFFFF'
  }),
  stroke: new ol.style.Stroke({
    color: 'black',
    width: 1
  }),
  text: new ol.style.Text({
    font: '18px Calibri,sans-serif',
    fill: new ol.style.Fill({
      color: 'black'
    }),
    stroke: new ol.style.Stroke({
      color: 'white',
      width: 0.5
    })
  })
});

function switchLabel(label) {

    if(label == 'a'){
        currentLabel = 'a';
    } else if (label == 'b') {
        currentLabel = 'b';
    } else {
        currentLabel = '';
    }

    var vectorLayer;
    var layersObject = map.getLayers();
    for (var i = 0; i < layersObject.array_.length; i++) {
        vectorLayer = layersObject.array_[i];
        break;
    }

    if(isFoo){
        vectorLayer.setStyle(function(feature) {
        style.getText().setText(feature.get(currentLabel));
        style.getFill().setColor(getColor(feature.get('pe_nr')));
        return style;
        });
    } else {
        vectorLayer.setStyle(function(feature) {
        style.getText().setText(feature.get(currentLabel));
        style.getFill().setColor(getColorByID(feature.get('id'), feature.get('pe_nr')));
        return style;
        });
    }
}


texts: Turning a single MultiPolygon into many small polygons?

texts:
I have a large GeoJSON file with a single very complex MultiPolygon, and I would like to turn it into an identical file with many small polygons.

How can I do this? I can use ogr2ogr or mapshaper or PostGIS, or any other command-line tool.

This is what my file looks like:


texts: Connecting ArcSDE geodatabse 10 from ArcSDE service Version 9.3

texts:
I have recently upgraded ArcSDE Geodatabse from 9.3 to 10.0.
The applications were using SDE connection through ArcSDE Windows service V 9.3

Now when I tried starting the same ArcSDE windows service on same instance and for same oracle database using following command - >>
sdemon -o start -i sde_pumaqa2 -p sde

This is giving following error -
sde_pumaqa2 service failed during initialization.
Please check event log or error log files.
Error starting sde_pumaqa2 service(0)
Could not start ArcSDE -- Check Network, $SDEHOME disk, DBMS settings and dbinit
.sde.

The event viewer shows following error -
The description for Event ID 110 from source sde_pumaqa2 cannot be found. Either the component that raises this event is not installed on your local computer or the installation is corrupted. You can install or repair the component on the local computer.

If the event originated on another computer, the display information had to be saved with the event.

The following information was included with the event:

sde_pumaqa2
Connection to RDBMS failed

the message resource is present but the message is not found in the string/message table

There was no change is any of the dbinit / driver file. They are checked and confirmed.
Please let me know if any step I am missing in this upgrade process.


texts:
I have previously drawn and added a SQLite line feature to a map with no problems. I get the attached message when I try to save it.  I can add and save as a shapefile but prefer to have only the one file. I was using Madeira before but using Coruna now in case that has any bearing on the problem.

As requested :

  Could not commit changes to layer derf

  Errors: ERROR: 1 feature(s) not added.

  Provider errors:
  SQLite error: unknown cause
      SQL: INSERT INTO "derf"("geometry",,"line") VALUES (GeomFromWKB(?, 4326),,?)

Have tried on Surface book with Madeira and it still works OK there.


texts:
I am facing a problem with the Layout API while implementing a plugin for QGIS 3. I created a layout template , which contains different elements like map, labels, legend and scale-bar.
I'd like to place my legend on the bottom of the page. But the size of the legend can change. So I wanted to place my legend using :
legend = QgsLayoutItemLegend(layout)
legend.setLinkedMap(map)
legend.setReferencePoint(QgsLayoutItem.LowerLeft)
legend.attemptMove(QgsLayoutPoint(page.pageSize().width()*0.79,page.pageSize().height()),useReferencePoint=True)

But this is not working and it gaves me the same result as without the setReferencePoint.
So does anyone know why is the "LowerLeft" not working ?
A solution to my problem would be to get the height of my legend and place it from the default UpperLeft position, but I wasn't able to find out how to get that. Any idea on how to get the size of an item ?


texts: Free high-resolution Orthophotography for Quebec

texts:
Is anyone aware of free high-resolution orthophotography for Quebec? Something 1m resolution or greater would be ideal.

I'm looking for imagery to use in a standard GIS as files (so not as WMS or other service, or via the web.)


texts:
I would like to be able to create a convex hull in ArcGIS Desktop 9.x, but I cannot find an appropriate tool.

How does one go about creating one?

I am interested in answers for all license levels: ArcView, ArcEditor and ArcInfo.


texts: Comparing string values in field calculator

texts:
I am trying to compare the values of 2 fields in my attribute table and update a 3rd field as a result. In the field calculator I have

CASE  when (("Address")  like  ("COUNTYNAME" )) then 'Match'
else 'Error'
End

so for the following toy data examples I would get the result Match,Error,Error
Address,COUNTYNAME

Smith, Smith County
Bob, Smith County
Jane, Bob County

I have also tried CASE  when (("Address")  IN  ("COUNTYNAME" )) then

Are there any other ways to compare text values within QGIS?


texts: Converting shapefile to CSV file using GeoTools in Java?

texts:
I am trying to convert my shape files to csv through a batch processing in Java using Geotools. The problem is I have more than 500 shape files (point features). The code is working fine for a small number of shp files. But when I am iterating the process over a large dataset (number of files), some of the files are not written and of 0 KB. I have no clue why this is happening.

I want all the files to be converted from shp to csv.

Here is my java code-

// shp2csv conversion
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.opengis.feature.simple.SimpleFeature;

public class GeoOpen
{

private static SimpleFeatureIterator simpleFeatureIterator;
static FileDataStore store ;
static String ID="";

public  static boolean openShapeFile(File srcfname) throws Exception

{
    String filename=srcfname.getAbsolutePath();
    File dataFile = new File(filename);
    dataFile.setReadOnly();
    store = FileDataStoreFinder.getDataStore(dataFile);
   // ShapefileDataStore store = new ShapefileDataStore(dataFile.toURL());

    SimpleFeatureSource source = store.getFeatureSource();
    SimpleFeatureCollection featureCollection = source.getFeatures();
    simpleFeatureIterator = featureCollection.features();
    return true;
}
public  static boolean iterate(File srcfname, String dest_path)
{
    File dest_file=srcfname;

    String dest_filename=GeoOpen.getFileNameWithoutExtension(dest_file)+".csv";

    try{
        BufferedWriter  bw=new BufferedWriter (new FileWriter(dest_path+dest_filename));

   bw.write("ID"+","+ "date"+","+"tstamp"+","+"X_prj"+","+"Y_prj"+","+"NEAR_FID"+","+"NEAR_DIST");
    while(simpleFeatureIterator.hasNext())
    {

        SimpleFeature f = simpleFeatureIterator.next();

     //  System.out.println(""+f.getID()+" ,"+f.getAttribute(1)+", "+f.getAttribute(2)); //ID, field 1: lat, field 2: long

   //writing in a CSV file
       bw.write("\n");

       int index = f.getID().lastIndexOf('.');
       ID=f.getID().substring(index+1); //get only the ID //by default ID comes with filename.ID format when getID() method called
      bw.write(Integer.parseInt(ID)+" ,"+f.getAttribute(1).toString()+", "+f.getAttribute(2)+","+f.getAttribute(3)+","+f.getAttribute(4)+","+f.getAttribute(5)+","+f.getAttribute(6));

    }

    }
    catch (IOException e)
    {
        e.printStackTrace();
    }
    catch(Exception ex){
        ex.printStackTrace();

    }
    finally{
        simpleFeatureIterator.close();
        store.dispose();
      }

    return true;
}

public static String getFileNameWithoutExtension(File f)
{
    String s="";
    int index = f.getName().lastIndexOf('.');
    //System.out.println (index);
  if (index>0&& index <= f.getName().length() - 2 )
  {
      s= f.getName().substring(0, index);
  }
  return s;
}
public static void main(String[] args) {
    try {

      String dest_path="C:\\ArcGIS test data\\commonfiles csv\\001\\nr_train\\";

      File src_dir= new File ("C:\\ArcGIS test data\\commonfiles prj\\001\\nearTrain");

      File [] file_nr = src_dir.listFiles();  //array containing files from nearRoad
        for (int i=0; i<file_nr.length; i++)  //near road
        {
            if (file_nr[i].getName().endsWith(".shp"))
            {
                GeoOpen.openShapeFile(file_nr[i]);
                GeoOpen.iterate(file_nr[i], dest_path );
            }
        }

    } catch (Exception e) {

        e.printStackTrace();
    }
}
}


texts:
I run gdaltransform like so:

 gdaltransform -s_srs nad27 -t_srs wgs84
 -104.5 39.5
 -104.5 39.5 0

 gdaltransform -s_srs wgs84 -t_srs nad83
 -104.5 39.5
 -104.5 39.5 0

 gdaltransform -s_srs nad27 -t_srs nad83
 -104.5 39.5
 -104.5 39.5000000009262 4.23276796936989e-005

In other words, if gdaltransform nad27 -> wgs84 -> nad83, the lat/lon do not change.

If, however, I transform nad27 -> nad83, my latitude is ever-so-slightly different.

Why, please? What is the correct process for transforming coordinates?

I have a legacy application that only saves the datum name (e.g. NAD27, WGS84, etc) and the lat/lon. In its processing, when trying to find the lat/lon in NAD83, which is the projection defined in a 1 sec FLT .HDR file, it calls (an old DLL implementation of) CS-Map cs_cnvrt, passing in nothing more than these datum strings, and it gets an output lat/lon that are different than the input lat/lon.

When updating this code, I have this crazy idea that I'll keep all coordinate internally in WGS84, and transform them as needed by the output plane. So my code reads the NAD27 coordinates, converts to WGS84 using GDAL. It later asks GDAL to transform WGS84 to the spatial reference system (SRS) in the FLT file (i.e. NAD83). The transform results are exactly the same as the input.

This makes sense to me, but I can't understand why the two methods above are so different.

Thanks.

-reilly.


texts: What is the proper name for the "Venn diagram" of a set of geometries?

texts:
Sometimes, you have a set of overlapping geometries, and you want to generate an output like this:

Image taken from this answer.

The resulting polygons have no overlaps; each region of space is covered by only one polygon. All boundaries are preserved; every combination of intersecting input polygons results in a separate polygon in the result. This is similar to a Venn diagram in that result is all the separate relations of the polygons.

What is the proper name for this operation or for the output? "Intersection" does not seem appropriate, as that might be only the region where all the inputs intersect. "Overlay" seems inappropriate, since that involves combining multiple sets of input (layers) and does not self intersect a single input layer.

I'm really at a loss at what to call this thing, but I feel confident it must have a formal name.


texts: Reading GeoJSON file and querying Nominatim?

texts:
I need to loop through a GeoJSON file, where for each "feature", find the city where the track started:

with open('input.geojson') as f:
    gj = geojson.load(f)

for track in gj['features']:
    #OK
    print track['geometry']['coordinates'][0]

    #How to query Nominatim to find the city that lat,lon belongs to?

    #Next step : edit 'name' in gj

with open('dump.geojson', 'w') as f:
   dump(gj, f, indent=2)

Would someone have some code handy?


texts:
I am attempting to create a binary cloud mask from the MOD35_L2 data products located at http://ladsweb.nascom.nasa.gov/data/ When opening the MOD35_L2 data in ENVI or the MRTSwath tool it reports 6 bands of data (which I believe is due to ENVI/MRTSwath reading the hdf data incorrectly) Is there a way to use ENVI or MRTSwath to create a MODIS cloud mask?

How can I process a binary cloud mask from the hdf files?


texts: ModelBuilder names not changing with batch operation?

texts:
I am using ModelBuilder for ArcMap.  I specify the names of my temp files like so:

%scratchworkspace2%\Ofst_B_%LetterId%_%Sea Level Rise Inches%.img

"LetterID" is a string parameter, "Sea Level Rise Inches" is a double, and "scratchworkspace2" is a folder.

When I run my model in batch mode these files are not stored with the correct name.  The first file is saved correctly I think, but the subsequent runs it saves the temp value as "G:\temp\rastercalc1".  It is not saving the temp file with the specified name.  I have to manually delete this file from the temp directory.

Also, the "Sea Level Rise Inches" value is not incrementing with the batch file, so it is overwritten with each batch run.

I am doing a lot with the % symbols and I may need some help with that.

But at an earlier stage of my model, I use a similar naming convention and the files do save correctly and are not overwritten:

%Working Directory%\Depth_Rast_%LetterId%_%Sea Level Rise Inches%.img

But this does NOT work without being overwritten:
    %Working Directory%\BC_%LetterId%_%Sea Level Rise Inches%.shp

Below is my model:


texts:
I have the task to provide my point observations data in the "WFS 2.0" (Web Feature Service 2.0) standard format.

Are there any free, easy-to-use, desktop GIS clients or online viewers for showing and querying my WFS 2.0 layers?

I have only found the "QGIS WFS 2.0 Client plugin"

Do any other popular free GIS viewers (Kosmo, MapWindow, ArcGIS explorer..) have support for showing and querying WFS 2.0 web feature service?


texts: Streaming polyline while adding records to related tables in ArcPad

texts:
I'm designing an ArcPad solution for our field staff; we collect point data and use lines to identify the amount of 'effort' (and to maintain accountability).

While collecting a 'streaming polyline' (Add verticies continuously from GPS) as our field staff's "track" (daily effort) we need to also collect point data along the way.

Every time you collect add a record to a point feature class in a related table in the poly-line streaming stops.  This forces our staff to manually turn the polyline streaming back on for every record they collect.  Human error is just about guaranteed over the course of editing 50 points in a day.

I know that event handling is built into ArcPad through ArcPad Studio, but I haven't had any luck triggering the polyline streaming with the onOK event of related table's form.

Any help putting together an event handler would be appreciated.  The help files aren't nearly as useful as the Python sample code associated with ArcMap on the desktop.  Here's a quick look at the direction I've been leaning. I'm not fully sure how to even compose a full even handling script for ArcPad.

  From help files:   'gpsstream' => Add Vertices continuously from GPS

  Form (onOk event):

  Application.ExecuteCommand("gpsstream")


texts: How to sequentially increment date field in ArcGIS

texts: I have created a date field in my point shapefile. I want to input date their. I am able to do so. Now I want to create a sequential increment with 6 second interval from previous point to next point. I have Google it but didn't find any solution. Please help me how can I do it with a huge data set.


texts:
I have create a application where we show world map using arcgis .Here user does lot of interaction with map he click on the map or we highlight map programmatically so i want to improve performance of map .

var WorldCountries = new ArcGISDynamicMapServiceLayer("urlr");
                map.addLayer(WorldCountries);
         map.graphics.enableMouseEvents();
                 featureLayer = new FeatureLayer("url",{
                    mode: esri.layers.FeatureLayer.MODE_SNAPSHOT,
                     autoGeneralize :true,
                     outFields: ["*"]
                  });

As below code shows i am using feature layer in snap shot mode so most of graphics are loaded in the client side,so i dont want to hit server again.So suppose in world map if i want to highlight some country i can use below code  so i dont want hit serve again.

dojo.forEach(featureLayer.graphics,function(graphic){
                                 if(graphic.attributes.CNTRY_NAME!=''  && alertNBulletinList.contains(graphic.attributes.CNTRY_NAME)){
                                     graphic.setSymbol(alertOfficeSymbol);
                                     map.graphics.add(graphic);
                                 }
                        });

   I have two question now :

1) is there a way i can cache loading of graphics , now suppose in landing page i load the map with all feature layer i go to other  page of my web application is there a way i can cache the loading of graphics since it is nothing but json object sent from server

2) Can i make use of feature collection to laod all the json object if so how do i query for some graphics


texts: Importing an image to QGIS with no coordinates, only a scale of the drawing

texts: Is it possible to import a drawining into QGIS without having any coordinates or CRS. The scale that I have is 1:1250.


texts: How to compute the smallest network that connects all points using PostGIS?

texts:
I have a set of postgis scripts which generates two tables - one of a set of points and the second a set of roads which surround them. All data is in the same projection and both outputs are stored in postgres 9.2 tables with postgis 2.1

The pgrouting topology of the road network has been created and the points table has a column containing the nearest road segment.

I'd like then to generate a subset of the road network which represents the smallest network that connects all points using something like a minimum spanning tree. The road network is undirected, and costs are simply the route length.

I can do this in QGIS/Grass using the v.net family of modules but ideally I'd like to keep this final step in SQL as well.

I've looked at the new apspWarshall postgis function but I'm at a loss as to how it can be encouraged to focus its energy on connecting the points and not the whole network.

This is the short script I've put together in an attempt to create a framework to solve this but I can't see where its possible to focus the function to start with a subset of the edges.

SELECT seq, id1 AS node, id2 AS edge, cost, the_geom
FROM   pgr_apspWarshall('SELECT gid AS id,
                                source,
                                target,
                                st_length(the_geom) AS cost
                         FROM   road_network
                        ',
                        false, false
                       ) AS tree
JOIN   road_network As roads
ON     tree.id2 = roads.gid

In single path shortest path problems the function asks for the start and end but apparently not in the all points problems. Equally in Grass the v.net.spanningtree and v.net.steiner expect a set of points and lines as a combined network to work with.

Does anyone have an suggestions for how to do this in PostGIS?


texts: Creating unique identifier for cadastral parcels in ArcGIS 10.5?

texts:
Creating unique identifier for cadastral parcels in ArcGIS 10.5 (ArcView)?

I'm adding selected parcels to a feature class and I would like to have a text prefix and then three or four digit number that increases as I add to the layer. The text prefix will be categorizing the type of parcel, the numbers will be for keeping a sequential count.

An old colleague gave me this script a few years ago, but no instruction.

Static recz As Long
Dim lStart As Long
Dim lInterval As Long
Dim IMap As String
Dim IMapID As String
'=================
'adjust start value and interval bellow
lStart = 500
lInterval = 1
IMap = "G505"

'=================
If (recz = 0) Then
  recz = lStart
End If
recz = recz + lInterval
IMapID = IMap & "." & recz


texts: Order of displaying labels within a layer

texts: How do I choose the order which labels are drawn, within a layer? I am using label classes but I don't want my points drawn in the order of the label class. My points are numbered and I want to ensure that the lowest numbered ones always sit above higher numbered ones.


texts: Converting a description of polygon with arcs to WKT

texts:
A document that I am trying to digitize contains descriptions of polygons (see below). The descriptions are converted to WKT for a backend. However, I get stuck at the "along (counter)clockwise arc" segments.
I have made several attempts to do this by having a CircularString that starts and ends at the specified points, and has the middle point at the radius, on the line from the center point of the circle, perpendicular to the line through the start and end points. However, this gave a very distorted shape.
For example, the bearing of the line between 53°06'18.47"N 006°24'23.79"E and 53°13'24.22"N 006°39'29.22"E is about 52 degrees. To get the perpendicular line (it is a circle, so some assumptions are made here), I set out a bearing of (52 deg - 90 deg) = 322 deg from the centerpoint of the circle. At 6.5NM on that line, I find the point 53°12'36.0"N 6°28'21.0"E. Drawn on a map, the segment looks like this compared to the intended shape:

I am currently stuck at how to approach this problem further, without converting the circle to a LineString first, and then intersecting with the start and end points. This is not my preferred approach because of the loss of precision.
53°02'27.95"N 006°16'17.31"E;
53°06'18.47"N 006°24'23.79"E;
along clockwise arc (radius 6.5 NM, centre 53°07'30.00"N 006°35'00.00"E) to
53°13'24.22"N 006°39'29.22"E;
53°17'13.14"N 006°47'39.82"E;
53°12'29.90"N 006°53'48.78"E;
53°08'40.58"N 006°45'36.79"E;
along clockwise arc (radius 6.5 NM, centre 53°07'30.00"N 006°35'00.00"E) to
53°01'35.60"N 006°30'32.01"E;
52°57'44.69"N 006°22'24.16"E;
to point of origin.

Note on this example is that the centerpoint for the circle is the same for both arc, this is not the case for every polygon.
The tools I have available are Python and QGIS. How should I approach this problem?


texts: Why does my 113th Congress Shapefile match up with the wrong districts across the country?

texts:
I downloaded a Shapefile from the Census on the 113th congress broken down by district nationally: http://www2.census.gov/geo/tiger/TIGERrd13_st/nation/tl_rd13_us_cd113.zip

Within the .dbf file, there were no good indicators of how the districts are normally laid out (Eg. PA-15); each state was assigned a number and then a congressional district, so I figured out what the state codes meant (based on other census shapefiles and how many districts were in each state and made a new column that I'm 99.9% sure is right based on the information given. The first four rows are what came with the shapefile and the last highlighted row is what I added:

Modified DBF.jpg

And here's what it looks like when I run a query for the districts that I've defined as Pennsylvania:

"Pennsylvania" GIS.jpg

Pretty sure that's not what Pennsylvania looks like! Can't figure it out -- because, again, I'm really sure that I defined each state's number as the proper actual state. Any ideas as to what might be wrong?

And sorry about all the links to pictures -- since this is my first post, I'm not allowed to post any images or more than two links yet.

Thanks!


texts: Networks plugin missing from QGIS 3.4.7?

texts:
I installed the Networks plugin in QGIS 3.4.7 successfully, but it does not appear in the Vector Tab.

Can anyone suggest another way to reverse the direction of a vector used in a network (graph)?


texts: Where can I find an administrative boundary shapefile for Monaco?

texts:
I am looking for an administrative boundary file (preferably in ESRI shapefile format) for the divisions of Monaco namely:

Condamine
Fontvieille
Jardin Exotique
La Rousse
Larvotto
Les Moneghetti
Monaco-Ville
Monte-Carlo

I have looked at a few websites like diva-gis, nationmaster, geofabrik, statsilk and eurostat but non of the files for Monaco are divided by the above areas, only the outline of Monaco itself. I need something akin to this image below...


texts:
I've got a data set that comes with coordinates, however they seem to have gotten scrambled somewhere along the way of data-collection, data-entering and data-export. For some of the places the actual coordinates are known. Now the goal is to use this to figure out how they got scrambled, in order to be able to reverse this.

Core question:

Which transformation has to be applied to data_lat and data_lon (below) to obtain coordinates in a known format.

The following table lists the approximate real coordinates (as decimal coordinates) and the scrambled coordinates I found in the data set. The latter  have an (or several) unknown format(s).

real_latiude,   real_longitde,  data_lat,   data_lon
13.30214,       -14.22023,      1318136,    141347
13.31396,       -14.21844,      1318268,    1412561
13.34097,       -14.09441,      1320262,    1405442
,               ,               1316060,    1408182
13.39157,       -13.88993,      568980,     1484519
13.43406,       -14.6541,       537144,     1485106
13.47477,       -14.08452,      599042,     1489640
13.48308,       -16.50089,      337657,     1491175
13.48558,       -15.96453,      395687,     1491387
13.48643,       -16.54477,      333733,     1491481
13.49866,       -14.82873,      518529,     1491920
13.59012,       -15.13228,      485738,     1502171
13.60742,       -15.13525,      485367,     1504185
13.61425,       -15.16174,      482322,     1504335
13.64409,       -15.10624,      488472,     1508243

1/3 of the solution:

For only the first three entries, I was able to figure out that the data was entered in UTC format and then all-special characters where stripped. I.e. it corresponds to the real coordinate, i.e.:

1318136,    141347 => 13°18'13.6"N  14°13'47"W => 13.303, -14.230
(this is close enough to the real_* coordinates)

Possible starting points

We have reasons to believe that the data was manually entered in Excel. So one possible thing could be that Excel interpreted the entered figures as dates somehow and save it's own weird date/time format. I was however so far not able to figure out how that could've happened.

Reward

I will set a bounty to this question (as many as I might have by then), as soon as I'm allowed to set a bounty.


texts: Accounting of strength of spatial autocorrelation between some grocery stores in distance 500 meters

texts: I'm analyzing sales one supermarket chain and I want to know how other neighborhood stores influence on turnover store. I have data about other stores and their location. As a result, I need to get the influence weights of neighboring stores on my stores. I am new to spatial econometrics but I read about Moran I and GWR. How can I calculate the value of communication strength between stores?


texts: Extent Intersection when WrapAround180 enabled

texts:
Allowing WrapAround180 allows map extents to go beyond the true Web Mercator extent. Once this happens I am not able to use map.extent.intersects(point). My work around was to

Normalize extents in lat/long
Convert point to lat/long
Split extents to two parts if crossing date line.
Test intersection with extent(s)

Is there an easier way?

Here's my code:

var extent = map.extent;
         var pointToTest = someGraphic;

     var normMin = wmUtils.xyToLngLat(extent.xmin, extent.ymin);

     var normMax = wmUtils.xyToLngLat(extent.xmax, extent.ymax);

     //test if minx is greated than maxx
     if (normMin[0]> normMax[0]) {

      var pointMin1 = new Point([normMin[0],normMin[1]],new SR({ wkid:4326 }));
      var pointMin2 = new Point([180,normMax[1]],new SR({ wkid:4326 }));
      var wmMin1 = wmUtils.geographicToWebMercator(pointMin1);
      var wmMin2 = wmUtils.geographicToWebMercator(pointMin2);

      var pointMax1 = new Point([-180,normMin[1]],new SR({ wkid:4326 }));
      var pointMax2 = new Point([normMax[0],normMax[1]],new SR({ wkid:4326 }));
      var wmMax1= wmUtils.geographicToWebMercator(pointMax1);
      var wmMax2= wmUtils.geographicToWebMercator(pointMax2);

      var normExtent1 = new Extent(wmMin1.x,wmMin1.y,wmMin2.x,wmMin2.y, new SR({ wkid:102100 }));
      var normExtent2 = new Extent(wmMax1.x,wmMax1.y,wmMax2.x,wmMax2.y, new SR({ wkid:102100 }));

    } else {

     var pointMin = new Point(normMin,new SR({ wkid:4326 }));
     var wmMin = wmUtils.geographicToWebMercator(pointMin);
     var pointMax = new Point(normMax,new SR({ wkid:4326 }));
     var wmMax = wmUtils.geographicToWebMercator(pointMax);
     var normalizedExtent = new Extent(wmMin.x,wmMin.y,wmMax.x,wmMax.y, new SR({ wkid:102100 }));

   }

   if (normMin[0]> normMax[0]) {
        if ((!normExtent1.intersects(pointToTest.geometry) && !normExtent2.intersects(pointToTest.geometry))) return;
      } else {
        if (!normalizedExtent.intersects(pointToTest.geometry)) return;
   }


texts: How to get values to create world file format?

texts:
i have researched regarding making world file format and tried to create one using c++ in windows, can somebody confirm if where i am getting the values from are correct? the literature regarding world file formats are quite vague. i don't use rotation so i assume value D and E are always zero..

Line 1: A: pixel size in the x-direction in map units/pixel.

p->GetDeviceCaps(HORZRES); // HORZRES is defined in wingdi.h as Horizontal width in
                           pixels

Line 2: D: rotation about y-axis // always 0 as i don't use rotation

Line 3: B: rotation about x-axis  // always 0 as i don't use rotation

Line 4: E: pixel size in the y-direction in map units, almost always negative

p->GetDeviceCaps(VERTRES); // VERTRES is defined in wingdi.h as Horizontal width in
                           pixels

Line 5: C: x-coordinate of the center of the upper left pixel
// how do you get this value?

Line 6: F: y-coordinate of the center of the upper left pixel
// how do you get this value?


texts: Distance in advanced digitizing not working properly?

texts:
I'm trying to remap a papermap of a festival with QGIS.
I tried to use the advanced digitalizing tool to set the right distances and angles between the point but noticed that the distance I put in isn't correctly set in the map.

To reproduce the problem do the following steps:

Load google satellite map
Go to coordinate X: 1041106.209718  Y: 7007221.582073 (this should be on a small airport in northern Germany)
Create a temporary layer with geometry set to polygon
Activate the advanced digitalizing tool and set the first point of the polygon at the given coordinate
Press D on the keyboard and enter any distance you want to (I used 25,611)
and set the second point and any third point so the polygon is valid
Use the measuring tool to check if the first segment has the desired length.
In my case it was permanently too short.

Does anyone know how to solve this problem?

I am using QGIS 3.8 on Linux Mint 19 Tara


texts:
How to automatically make the title of your print equal to the name of the printed layer? Tried by using the next piece of code, but it doesn't work (QGIS 3.X):

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

@qgsfunction(args='auto', group='Custom')
def func(feature, parent):
    root = QgsProject.instance().layerTreeRoot()
    layer = qgis.utils.iface.activeLayer()
    lyr = root.findLayer(QgsProject.instance().mapLayersByName(layer.name())[0].id())
    group = lyr.parent().name()
    if group != '':
        return lyr.parent().name()
    else:
        return 'Group not found'


texts: ST_Intersects PostGIS returns empty table

texts:
Here is a link to download the two layers I am working with (as geopackages).

I have two tables in a PostgreSQL database, 'WataugaResidential' and 'zone1shape'. I'm trying to use ST_Intersects function from PostGIS to see where they overlap. Using QGIS, if I run the Vector > Research Tools > Select by Location I get the following result:

Trying to do the same in PostgreSQL, I get an empty table returned as my result. Here is the SQL code I'm using:

SELECT res.*
FROM public."WataugaResidential" as res, public."zone1shape" as z1
WHERE ST_Intersects(res.geom, z1.geom);

Pretty simple stuff. But my result is an empty table. All the columns are there, but no data. Result:

Why is nothing returned?


texts: Bad module reference in GDAL -- QGIS 3.10.1

texts:
Trying to use GDAL with QGIS 3.10.1, I encounter the following error:
GDAL command:
python3 -m gdal_polygonize contextualWMSLegend=0&crs=EPSG:4326&dpiMode=7&featureCount=10&format=image/jpeg&layers=25&styles&url=https://maps.vcgi.vermont.gov/arcgis/services/VCGI_services/VCGI_BASEMAP_SP_v2/MapServer/WMSServer "C:/Users/Dan Grossman/AppData/Local/Temp/processing_53de98b33ff64a33948c3e46ee96da46/24da5c228f344e2a8c113d3c7a896890/OUTPUT.gpkg" -b 1 -f "GPKG" OUTPUT DN
GDAL command output:
C:\PROGRA~2\QGIS3~1.10\bin\python3.exe: No module named gdal_polygonize

'crs' is not recognized as an internal or external command,
operable program or batch file.
[etc.]

gdal_polygonize.py is in the QGIS 3.10 > apps > python37 > scripts directory, so I have to assume that the module did not get include(ed). Other GDAL algorithms have the same problem. Beyond that, I'm stumped.


texts: How to spatially join OSM point to road data?

texts:
Using QGIS 2.16.1 , using shape layers, with the same CRS, I am having two  line and a point layer.  The point layer was added by the plugin : locate points along lines" The third layer is the downloaded layer "roads" (shape-OSM) from geofabrik. I want to add the streetnames from the roads layer to the points layer. The points are on te line.....  If possible without extra steps like buffers..

Which spatial join should I use? And, do I use intersect/touches/crosses etc?


texts:
While importing TIFF image to GRASS GIS, I get an error:

Location PROJ_INFO is:
name: Lambert Conformal Conic
proj: lcc
datum: nad83
a: 6378137.0
es: 0.006694380022900787
lat_1: 36.16666666666666
lat_2: 34.33333333333334
lat_0: 33.75
lon_0: -79
x_0: 609601.22
y_0: 0
no_defs: defined

Dataset PROJ_INFO is:
name: Lat/Lon
proj: ll
datum: wgs84
ellps: wgs84
no_defs: defined

ERROR:
You can use the -o flag to r.in.gdal to override this check and use the location definition for the dataset.
Consider generating a new location from the input dataset using the 'location' parameter.
Provide rolution to remove this error -o in GRASS GIS


texts: Getting ArcGIS online web map to update images coming from randomly generated link?

texts:
I am creating a web map that shows images from a live traffic feed. The feed, which is coming from the Texas Department of Transportation's website, updates every few minutes. When the updates happen on TxDot's website, a link is randomly generated.

That means when I copy and paste the link from the website into my CSV, I am pulling one of those randomly generated links. The links will only show the picture from when the link was pulled and will not update with the most current picture from TxDot's website.

Is there a way, either within the CSV file I am building to create the web map or within ArcGIS online, to get the most up-to-date link without having to manually input it? Each image I am trying to capture is uniquely named within the TxDot website. That unique name is used in my CSV.

I want to clarify that I have no problem getting my CSV into ArcGIS Online. I have a complete web map that looks great. My issue is that my pictures will not update because, essentially, the links expire. Is there a way to build a script that can access a website and update the link based on image name?


texts: Creating custom tool using Python script exported from ModelBuilder?

texts:
I created some ModelBuilder for some customization to reduce time in our routine process.

I exported my ModelBuilder to python script by availability of export options.

It doesn't work if i bring to tool bar.

any one suggest me clearly how to access the script in tool bar with some standard clear steps to follow?

as well as whether i need to do any modification on my exported python script?


texts: QGIS 'Voronoi polygons' freezes with 500k points

texts:
I've been using the 'Vector geometry--> Voronoi ploygons'function in QGIS to create voronoi cells for 200,000 points with success. However, recently I am working with a bigger point set (500,000 points randomly generated into a polygon repreesenting 1000km2) and the function consistently freezes without error at 99%.

I have tried outputting both to a temporary file (memory) and to a package file, but both still freeze at 99%.

Each point has coordinates and several attributes, so I thought perhaps the file was just too large for QGIS. I have tried running again with only the coordinates, but it still hangs at 99%.

Can anyone advise?


texts: Cartesian co-ordinate system with set 0,0 point?

texts:
This is for visualization in construction project planning:

I would like to set a 0,0 point and angle relative to North for an X,Y (and possibly Z) Cartesian co-ordinate system that will then show the cursor position in feet (or meters) relative to the 0,0 point.

It would behave like the co-ordinate window showing Lat./Long. but in feet instead, using the projection system of the layer, or project.

A specific example: I would like to set the corner of a barn as 0,0 and then know where the cursor is in X,Y feet (aligned to the property's road frontage) from that point as I move it around the map.


texts: How do I configure GeoWebCache to read a TMS (GoogleMaps or OSM)

texts:
Possible Duplicate:
  geoserver support to proxy other tiling services?

Probably a newbe question, with a simple answer, but I can't find enough detail in the GWC documentation.

I would like to set up GWC so that it reads a TMS (At the moment, just GoogleMap or OSM) and then makes this available as a WMS layer through GWC.

Is this possible, and if so, where do I configure it? I have been looking through the documentation on how to edit and add layers in geowebcache.xml, but I can't see the extact syntax for adding a TMS.

Any help would be greatly appreciated!

Thanks


texts:
With QGIS 3.8 a new Tool has been introduced: Overlap Analysis. It says that it "calculates the area and percentage cover by which features from an input layer are overlapped by features from a selection of overlay layers".

Now I have a feature layer including spatial data of operating areas of my company (which unfortunately are not congruent with administrative regions), which I would like to intersect with a layer of administrative areas of Germany (landkreise) and some other information which I only have for these administrative regions. As a result I need a new layer which shows me to what degree my company's operating areas are overlapped by Germany's administrative regions (landkreise).

So I chose one operating area and the overlapping administrative regions on the map. As Input Layer I chose my companys operating areas and as overlay I chose the administrative regions. Once I run the tool, the result is not as I imagined it. Most of the time I only get one feature in the resulting layer, but it doesn't show me to what degree the administrative regions overlap the operating areas.

Anyone got any experience with the new tool overlap analysis?


texts: Manually Adding a Decimal Field to the Attributes Table

texts:
I am fairly new to QGIS and search as I might in the manual, training manuals and widely on the net I can't find an in-depth explanation of how to set up a Decimal Number (real) field in my attribute table properly.

I have a number of vector layers, each with a polygon on them.  I need to add a decimal field to manually insert a rate per 1,000 of the population.  The number is more than likely going to be between n.n and nn.n

I have tried all sorts of different width and precision figures, but either the field is too small to take the number, or it is large enough, but I can't get a decimal point in there.  Does this mean that the decimal number is less than one?

When I have entered some information into the cell and pressed Enter the cell always reverts to NULL.

Please help!

Thanks

Andy


texts: Can QGIS be used behind paywall or barrier page?

texts:
The company I work for is a pay-for-subscription news website.

Is it legal to use QGIS to host maps that would be only available on our site to paying subscribers?

I'm having trouble figuring out the GNU license.


texts: Serving GroundOverlay image data from a KML file via WMS

texts:
I've got a KML file with a collection of sparsely distributed (and small) GroundOverlays, with each overlay containing a distinct LatLonQuad and Icon.

I would like to serve the this image data via a WMS (ideally GeoServer). The Icons (images) in the KML are aerial photos of real-world objects, so it is imperative that they scale correctly according to the zoom level of the view.

I was thinking of writing a script to convert the KML file to one or more GeoTIFF files or something similar (as described here: Export georeferenced image overlay from KML?), but it feels like that is not the most efficient approach, particularly considering that the GroundOverlays are so sparsely populated.

Another option I was looking at was trying to just create a Shapefile from the KML file and use SLD with a dynamic symbolizer (as described here: https://stackoverflow.com/questions/2589181/display-image-over-map-layer-in-geoserver) to avoid creating a raster for such sparse data, but I suspect that the icon images will be tiled instead of stretched at different zoom levels, which is not the intention.

Is there a quick and easy way to serve this kind of data via a WMS?


texts: How to draw on an ESRI.PrototypeLab.Surface?

texts:
i'm trying to draw on the ESRI map enabled to Microsoft Surface (ESRI Map -> Surface) but i can't find the way to do so. The problem is that the Draw method (which come from the ArcGIS API for WPF, the original one) needs a map ESRI and with this Surface project, this type of map is different (necessary to run on Surface). I thought that using Reflector, i could open the original Draw method and try to modify it to make it work with the "Surface Map" but it's way too complex for me ... Does someone already tried to do that kind of thing ? I'm really stuck with it! If it's not clear enough, just let me know. Thanks in advance.

edit :

private void SurfaceButton_ContactEnter(object sender, ContactEventArgs e)
    {
        _graphicsLayer = surfaceMap.Map.Layers["MyGraphicLayer"] as GraphicsLayer;

        var pointConverter = new PointCollectionConverter();
        var pointCollection1 = pointConverter.ConvertFromString(_data) as ESRI.ArcGIS.Client.Geometry.PointCollection;
        var polygon1 = new ESRI.ArcGIS.Client.Geometry.Polygon();
        polygon1.Rings.Add(pointCollection1);

        var graphic = new Graphic()
        {
            Geometry = mercator.FromGeographic(polygon1),
            Symbol = LayoutRoot.Resources["DefaultFillSymbol"] as Symbol
        };

        if (_graphicsLayer != null)
            _graphicsLayer.Graphics.Add(graphic);
    }

Here is the code of my button, which should add a polygon into the map. But nothing appends.
On my graphicLayer : Opacity = 1.0, Visibility = true. And the DefaultFillSymbol has Red properties.

edit : Maybe something wrong with my layer ?

GraphicsLayer gl = new GraphicsLayer()
                               {
                                   ID = "MyGraphicLayer",
                                   Opacity = 1,
                                   Visible = true
                               };

        surfaceMap.Map.Layers.Add(gl);

        if (!gl.IsInitialized)
            gl.Initialize();


texts:
I am having issues with a data layer that is not lining up with my other data layers!  What do I do?

I have all the parcel polygons, parcel pts. Jurisdiction polygons from Ashtabula county that all line up (I received these data layers all from the same source...from Ashtabula county).  Then I bring in a data layer for street centerlines from OGRIP (where I have obtained all my other street centerline data from that works great!)

How do I get these to work?  The extents on the street centerline data are completely different from the other data layers.

The projected coordinate systems are the same in all the layers?


texts: Calculating limits of visible horizon of location in QGIS?

texts:
I am an amateur at GIS. I installed QGIS Las Palmas version in Windows, I searched YouTube and also in Google but I didn't find this question.

There is a website that one can input the coordinates of a location and its elevation above ground level, then the website calculates the viewshed (cloak visibility) and the horizon limit as a kml polyline.

Having the DEM (.hgt file) for the desired region, how can I calculate the real horizon limit (as .kml polyline/polygon), not the viewshed, taking into account terrain in QGIS, just the limit of visible horizon (that it would be theoretically a circular polyline around the location without accounting for terrain).

Such a thing:

I would have such a horizon with QGIS, because the rendered horizon from the the cited website is not for each compass bearing, I am looking for a more accurate horizon line.

My goal is to compute at each bearing, the minimum angular elevation of sun or other objects in order to be seen from the location of observer.


texts: Importing shapefile into PostgreSQL failing

texts:
I have installed PostgreSQL 12 on my Windows 10 OS. I installed the PostGIS extension as well. I would like to import a shapefile and am
using shp2pgsql-gui to do so. The import is failing and there is no explanation why. Any ideas?

I should also add that my bigger objective is simply to load US census block shapefile data into a PostgreSQL database that is hosted by AWS RDS. I figured the easiest way to do this is to download shapefile data to a db instance on my computer and then dump this database into the RDS database. If there is an easier way, I'd accept that as well.


texts: Getting error in applyEdits function for ArcGIS

texts:
I have been trying to add an entry to a featureLayer and in want to get a callback from this so I can update another field in it. below is my code

clkListenerAddFeature = dojo.connect(map, "onClick", function (evt)
{
    //alert("inside the click function");
    var point = evt.mapPoint;

    var attributes = {};
    attributes["1"] = "CLOSED";
    attributes["2"] = "DENIED";
    attributes["3"] = "0";
    attributes["4"] = "1";
    attributes["5"] = userLoginId;

    var graphic = new esri.Graphic(point, null, attributes);
    featureLayerEdit.applyEdits([graphic], null, null);
    dojo.connect(featureLayerEdit, "onEditsComplete", addResultsComplete);

    if (clkListenerAddFeature) dojo.disconnect(clkListenerAddFeature);
    refresh_map();
    reset_map(map);

});

function addResultsComplete(addResults)
{
    if (addResults.length > 0) {
        var graphic = featureLayerEdit.graphics[featureLayerEdit.graphics.length - 1];

        graphic.attributes["6"] = graphic.attributes["OBJECTID"];
        featureLayerEdit.applyEdits(null, [graphic], null);
    }

}

but I keep hitting this error TypeError: Unable to get value of the property 'getCellsInExtent': object is null or undefined

anyone has any ideas?


texts: Rendering a QGIS atlas from PyQGIS from outside of the main application

texts:
I'm struggling to get QGIS to render an atlas from PyQGIS outside of the main application. Version is a recent git master checkout. The coverage layer is loaded, multiple images are rendered but none of the composer label expressions based on the coverage layer fields evaluate (I get syntax errors not seen in the main interface) and the rendered images do not have the atlas-driven map. The non-atlas driven overview map is rendered but does not show the atlas coverage overview.

Atlas export works fine from the GUI application using the same project.

I've tried a few permutations of the code below but I'm not making much progress. Any clues as to what I'm missing here?

Sample syntax error:

Found expression: UTMMajorLabel(round(Y_min(map_get( item_variables( 'main'),'map_extent')),0))
../src/core/expression/qgsexpression.cpp: 441: (replaceExpressionText) [0ms] Expression parser error: syntax error, unexpected '(', expecting $end

Python script:

import sys
from qgis.core import (
    QgsProject, QgsComposition, QgsApplication, QgsProviderRegistry, QgsLayoutManager, QgsMapLayerStore, QgsAtlasComposition, QgsComposerMap)
from qgis.gui import QgsMapCanvas, QgsLayerTreeMapCanvasBridge

##Set these
#export QGIS_PREFIX_PATH=/usr/local/Cellar/qgis3-dev/2.99/QGIS.app/Contents/MacOS
#export PYTHONPATH=/usr/local/Cellar/qgis3-dev/2.99/lib/python3.6/site-packages/

gui_flag = False
QgsApplication.setPrefixPath("/usr/local/Cellar/qgis3-dev/2.99/QGIS.app/Contents/MacOS", True)
app = QgsApplication([], gui_flag)
app.initQgis()

project_path = 'blahblahblah.qgs'

canvas = QgsMapCanvas()
project = QgsProject.instance()
project.read(project_path)
bridge = QgsLayerTreeMapCanvasBridge(
  QgsProject.instance().layerTreeRoot(), canvas)
bridge.setCanvasLayers()

layoutManager=project.layoutManager()
comps=layoutManager.compositions()
composition=layoutManager.compositionByName('Print1')

atlasMap = composition.getComposerMapById(7)
print("Atlas Composer Map: {}".format(atlasMap.displayName()))

atlas = QgsAtlasComposition(composition)
coverage = project.layerStore().mapLayersByName("coverage")[0]
atlas.setCoverageLayer(coverage)

atlas.setEnabled(True)
atlas.setHideCoverage(False)

atlasMap.setFrameEnabled( True);
composition.setAtlasMode(2)

atlasMap.setAtlasDriven(True)
atlasMap.setAtlasScalingMode( QgsComposerMap.Auto )

atlas.beginRender()
print('{} Atlas Features'.format(atlas.numFeatures()))

for i in range(0, 3):
     print('page: {} {}'.format(i,atlas.nameForPage(i)))
     atlas.prepareForFeature( i )
     atlas.
     image = composition.printPageAsRaster(0)
     fname=('{}.jpg'.format(i))
     image.save(fname)
atlas.endRender()

QgsProject.instance().clear()
QgsApplication.exitQgis()


texts: Create new raster from both overlapping and non-overlapping values of two other rasters in R

texts:
I am having trouble figuring out the best way to create a new raster from values in two separate rasters.  Some solutions involve creating stacks or bricks, others do not.

Here is an example:

First raster is "propSummPCPmod4"

Determine frequency of values

freq(propSummPCPmod4)

 value    count

[1,]     1  1417084

[2,]     3  9103585

[3,]    NA 14535331

Second raster is "propWintPCPmod4"

Determine frequency of values

freq(propWintPCPmod4)

 value    count

[1,]     2   241512

[2,]     3 10279157

[3,]    NA 14535331

create a brick and plot both rasters side by side

brick1 <- brick(propSummPCPmod4, propWintPCPmod4)
plot(brick1, col=terrain.colors(255))

The plot shows that values of 1 from "propSummPCPmod4" [layer 1] and values of 2 from "propWintPCPmod4" [layer 2] do not overlap.

I want to a create a new raster and to employ the "1" values from propSummPCPmod4, the "2" values from propWintPCPmod4, and to have all other values represented by "3" or NA.  (NA values are outside the boundaries of the original source rasters)

I tried merge, cvr, overlay, and mask but couldn't get a satisfactory result.


texts: Using live feeds from ArcGIS Online for analysis?

texts:
I would like to import into arcGIS Desktop a livefeed layer and use it to perform some analysis.

This is the link to web page: http://tmservices1.esri.com/arcgis/rest/services/LiveFeeds/Shakemap/MapServer

If I download the layer into ArcGIS for Desktop, this is not editable.

Does anyone know a solution to this?


texts:
I have got the forest cover data of a region in africa (raster layer) and some data collected inside this region by a survey (.csv). In the .csv file I have got x,y coordinates and 3 other columns with information.

The coordinates are in formated in WGS84 projection as my raster file. Now i need to add 1 point for every coordinate into the raster layer.

I guess my problem is mainly that I have no allocation to a CRS in my .csv data.

How can I solve that or is that not the problem?


texts: Consolidating many-to-one overlapping polygon feature values into another polygon attribute using ArcGIS for Desktop?

texts:
I have a situation similar to the one here:
Using Spatial join for many to one relates in ArcGIS for Desktop?

however the desired end product is a dissolved/aggregated feature that collects up all the values it was joined to.

For example if I have PolygonA which is overlapped by Value1 and Value2 I want to produce a polygon feature class with the same geometry as PolygonA, but with an attribute value that has "Value1, Value2".

I have found some variations on answers for this, but hoping for something simple-ish that doesn't require a bunch of python (not a coder).  Hoping for a geoprocess that could handle the aggregation after a spatial join, does anything like that exist, or is there a technique to accomplish this?


texts: ERROR 999999 when mosaicing GeoTIFFs in ArcGIS

texts:
I have been trying to use the mosaic tool in ArcMAP to merge several GeoTIFFs together, using BLEND as Mosaic Operator and LAST as Mosaic Colormap Mode. I placed all files that I want to merge in the same folder, and also copied the first one as a target file.

That way, I had 8 folders that needed to be mosaiced. For the first two I had little issues, but with all other the program kept giving me an error message.

I tried various things, nothing worked.

How can I resolve this error?


texts: Error open KMZ file from Oregon 650t Garmin

texts:
This is similar to this Previous post however; it does not mention about dpi or resolution to save a good quality image with very limited kmz files to be use for Garmin?

What is the dpi when do you save the kmz file from Google Earth Pro into Oregon 650t Garmin ? Is there a standard file number  to save a kmz  file with a specific number ?

I had a Non- user  set up his ArcMap 10.1 and export user's map to jpeg and gave them to me to geoference them in GEP .

My question Do I need to know the Map scale from ArcMap exported to jpeg file ? or regarding of what the user map scale set at and export jpeg at a higher resolution to preserve a good quality image ?

Garmin will not open kmz files greater than 1 MB ? Am I wrong ? or do I need to resize the pixel size ?

My Garmin GPS does not have usb chip in, so I am thinking I might need that for the CustomMaps folder ?


texts: Zonal Statistics giving ERROR 000824: The tool is not licensed

texts:
I have a script and when I run the ZonalStatisticsAsTable I get this error:

File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\geoprocessing_base.py", line 510, in
return lambda *args: val(gp_fixargs(args, True))
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000824: The tool is not licensed.
Failed to execute (ZonalStatisticsAsTable).

My code is:
flist = []

    for file in files:
        if os.path.isfile(os.path.join(root,file)) and file.endswith('.shp'):
            flist.append(os.path.join(root,file))```

    arcpy.env.workspace =r"Z:\EMSV-057_AOIs_Reshape\5_Process\1_TaskA\ModelBuilder\Pruebas_MV\EU_AOI\CAPAS"
    for z in flist:
       shape=z

       erase=arcpy.Erase_analysis(shape,"C:\\OSM_coastlines\\water_polygons.shp", "Tierra", cluster_tolerance="")

       tZONA=arcpy.gp.ZonalStatisticsAsTable_sa(erase, "FID", raster, "Tabla", "NODATA", "ALL")


texts: How to link spatial and non spatial table in geoserver

texts: I have a vector data published in GeoServer. I have other two non-spatial table which can be linked to the vector by a common field name. Now I would like to know how this can be done using GeoServer?I had read something regarding views in GeoServer.Will the views helps in linking this?


texts: Finding nearest point to a polygon using QGIS

texts:
I am using QGIS.

I would like to get the closest points to each building I have already selected. I have layer for points (low tension poles) as well as layer for buildings, I selected some buildings which are quite many and wish to get the closest point (pole) to each building I have selected.

Is there any plugin or faster process I can use to solve the issue.


texts: Plot only the borders of the states excluding the national boundary

texts: Is there a way I can just plot the boundaries if the states excluding the national boundary of the US? I have both nation as well as state shape files from the census data. I would preferably like to use simple features (sf package) in R.


texts: Using Python in the Field Calculator, how do I refer to a field in another shapefile?

texts:
I want to compare the changes in geometry between two shapefiles (of election districts in NYC. One is from 2001 and the other 2005). I am attempting to calculate a new field in the 2005 shapefile that shows the change in AREA from the 2001 shapefile to the 2005 shapefile, using the districts common IDs.

What I can't seem to figure out is how to refer to a field name in ANOTHER shapefile in my field calculator definition.

Here's what I have so far:

Pre-Logic Script Code:

def areaChg(shpOrig, areaOrig, shpNew, areaNew):
    origRows = arcpy.SearchCursor(shpOrig)
    for origRow in origRows:
        newRows = arcpy.SearchCursor(shpNew)
        for newRow in newRows:
            if origRow.ED == newRow.ED:
                return newRow.ED - origRow.ED
            else: pass

AreaChg =

areaChg("nyed_01bav", areaOrig, "nyed_05", !Area!)

It's the "areaOrig" in that last part of code, which is the place-holder for the field in the other shapefile, which I'm having a hard time with. How do I point to a field in another shapefile in field calculator? For the field in the shapefile from within which I'm doing the calculation, the correct format is "!Area!", but what about for a field from another shapefile?

Thanks!


texts: Google maps api, several moving markers using a csv file containing the coordinates

texts:
I collected 5 sets of coordinates of moving vehicles and I want to play the scene using googlemaps API showing those 5 vehicles moving.
The sets are saved in a postgresql table.

is it possible ? if yes how ?

what I found for the moment is only one vehicle.


texts: Extracting info on polygons from underlaying raster layers

texts:
In ArcGIS 10.1, I have two layers:
1) raster of Topographic Wetness Index
2) vector layer of polygons representing soil types

I am wondering how to accomplish what I am after, that is to get information for polygons from the underling raster. In other words, I would like to extract the mean TWI for each polygon in order to calculate  what is (on average) the TWI for each soil type.


texts: Increase elevation for a portion of the raster covered by a shapefile

texts: I have a raster of California. I have the shapefile over CA. I would like to keep the bits hanging off into OR, NV, and AZ. What I'd like to do is to raise CA up about 200 feet from where it is, making it "pop out" from the surrounding states. There has got to be a way to use the Serval toolbar to select the raster portion covered by the shapefile and then use the raster calculator to add x-feet to that selection.


texts: Getting distance between two points over linestring using PostGIS?

texts:
I have a linestring:
LINESTRING(-1.3326397 50.9174932,-1.3319842 50.9166939)

and I select two points that lay on that linestring.

How do I get the distance between these two points (in meters)?


texts: CQL Filter for geoserver geojson service with fields with space

texts:
I have trouble with CQL filter in geojson.

The CQL filter that I want:

"Sector Catastral" = 'VERONA'

The URL that compose for this purpose:

http://localhost:3000/geoserver/SIGUD/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=SIGUD:vista_lotes&outputFormat=application%2Fjson&CQL_FILTER=%22Sector%20Catastral%22%20=%20%27VERONA%27

Response:

This XML file does not appear to have any style information associated with it. The document tree is shown below.
<ServiceExceptionReport version="1.2.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd"><ServiceException>
  Could not parse CQL filter list.
Encountered "\"" at line 1, column 1.
Was expecting one of:
    "not" ...
    "id" ...
    "in" ...
    <IDENTIFIER> ...
    <DATE_TIME> ...
    "(" ...
    "[" ...
    <DATE> ...
    "-" ...
    <INTEGER_LITERAL> ...
    <FLOATING_LITERAL> ...
    <STRING_LITERAL> ...
    "true" ...
    "false" ...
    "point" ...
    "linestring" ...
    "polygon" ...
    "multipoint" ...
    "multilinestring" ...
    "multipolygon" ...
    "geometrycollection" ...
    "envelope" ...
    "include" ...
    "exclude" ...
     Parsing : "Sector Catastral" = 'VERONA'.
</ServiceException></ServiceExceptionReport>

The CQL filter that actually works:

"Localidad"='CHAPINERO'

The URL that compose for this purpose:

http://localhost:3000/geoserver/SIGUD/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=SIGUD:vista_lotes&outputFormat=application%2Fjson&CQL_FILTER=%22Localidad%22=%27CHAPINERO%27

Response normally with filtered features.

I want to filter by field with spaces.

My geoserver version is 2.7.6

With OGC Filters i'm getting the same error by the space in field name.


texts: How to computational efficiently check if there are any overlapping polygons in a group of polygons?

texts:
I need to do some processing with groups of polygons. If any of the polygons in the group overlaps other polygon the result is unpredictable and don't throw any error. So, I need to detect if there is any overlapping polygon inside the group and which polygons are overlapping.

The groups are made of about 15 polygons, so checking each polygon against the others would be costly. I wonder if there is an efficient way to do that.

The code is in Python and can use either the OGR bindings or Shapely.

A theoretical answer is also appreciated.


texts: Using FME to write to variable output files?

texts:
In FME (2016.1), I am combining data from 10 shape files to one big layer, and want to write those to a DGN file. But only for a variable area name. There are more area names in the combined shape files, and per area I want to create a DGN with that area name as filename.

Can this be done in FME?


texts: Open raster referenced in attribute column of vector layer

texts: I have a polygon layer. One of the columns holds an absolute path to a (georeferenced) raster. Is there a way to (programmatically) open this raster in QGIS and add it to the map canvas?


texts: What is hierarchy (object model) of features in ArcGIS?

texts:
I'm new to GIS.

My understanding of shapefiles is that they are (contain?) featureclasses, that contain features (points, polylines, polygons). Is there a clear hierarchy for these files. Would the next level after features be individual points? Are polygons and polylines merely point collections? What implication does this have for working with the data - e.g. can all geometries be worked with at the "point" level?

Is there a resource that would help me to understand the data structure and/or object model of these files? What I am essentially looking for is an explanation as to what exactly shapefiles are and what I should know when working with them, especially when using ArcObjects. In addition to this - how does this relate to other software like QGIS? Do they operate in the same way?


texts: Help Setting up GeoServer Proxy with OpenLayers

texts:
I have a GeoServer WMS in Apache-Tomcat at http://localhost:8080 and a mod_wsgi app in standard Apache at http://localhost:80/myapp.

My OpenLayers map is being hosted statically by the Apache-Tomcat Geoserver app at http://localhost:8080/geoserver/www/OL.html. The JavaScript in this map has to call the mod_wsgi app for non geographic data.

I've tried setting up a proxy in GeoServer the following ways

http://localhost:80/mypp

http://localhost

localhost

None of these work. After making a request some of my response object attributes look like this which means it's not really working. GeoServer logs show no errors

priority "NORMAL"
readyState 4
responseText ""
responseXM Lnull
status 0
statusText ""

Can't say I'm too much of an expert at Proxies. But it seems to me that I wouldn't need to set up a OpenLayers proxy (as mentioned here) if GeoServer can proxy for me.

Ideas, hints, advice?


texts:
I am trying to select line features from an OSM road shapefile layer (Layer 1) based on the features of a line layer (Layer 2). Layer 2 is selected road segments from Layer 1 produced using a routing plugin in QGIS

What I want to do now is to select the road segments from layer 2 to layer 1. I am using the "Select by Location" tool on ArcMap using "intersect" as a spatial selection method but it selects segments that I do not want to be included.

I tried to use other spatial selection methods but to no avail- they do not select exactly the segments I want. The picture below shows "are within the source layer feature".


texts:
In the picture you can see, the red lines i added manually. These are the ends of the lines. To distinguish that there is the beginning of a new river segment, i would like to have something like these read breaks at the end or between each line segments automatically.
Is there any posibility to do this?


texts: Marker differently depending on the tool with Leaflet.Editable

texts:
I created a toolbar to add different markers with Leaflet.Editable (http://leaflet.github.io/Leaflet.Editable/example/index.html)
Is there a way to change the marker according to the tool used to add a marker?
How can I pass the attributes to be a L.Control.extend to a layer designed by Leaflet.Editable?

Any suggestion is welcome!


texts: Providing multiple shapefiles belonging together in single file with predefined symbology in ArcGIS Desktop?

texts: There are multiple shapefiles which belong together and should be provided in a single file with predefined symbology (see figure). Is there an option to do so? Currently, the shapefiles are bundled in a geodatabase file (red box in figure). However, it seems to be impossible to define a symbology in a geodatabase file. Another option is to produce a layer file (.lyr) which then holds a link between the shapefiles and the desired symbology. The drawback then is, that many files (shapefiles, layer file) have to be kept together.


texts: What does the "scale" setting when plotting a NDVI or EVI time series?

texts:
I adapted the script from the guide for ui.Chart.image.series to my region of interest and time period:
// Import the example feature collection and subset the forest feature.
var bm = ee.FeatureCollection("users/fulano_letal/ROI_BM_Mayor100ha_Veg1aria");

// Load MODIS vegetation indices data and subset a time period of images.
var vegIndices = ee.ImageCollection('MODIS/006/MOD13Q1')
                     .filter(ee.Filter.date('2000-01-01', '2017-12-31'))
                     .select(['NDVI', 'EVI']);

// Define the chart and print it to the console.
var chart =
    ui.Chart.image
        .series({
          imageCollection: vegIndices,
          region: bm,
          reducer: ee.Reducer.mean(),
          scale: 50000,
          xProperty: 'system:time_start'
        })
        .setSeriesNames(['EVI', 'NDVI'])
        .setOptions({
          title: 'Average Vegetation Index Value by Date for BM',
          hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
          vAxis: {
            title: 'Vegetation index (x1e4)',
            titleTextStyle: {italic: false, bold: true}
          },
          lineWidth: 5,
          colors: ['e37d05', '1d6b99'],
          curveType: 'function'
        });
print(chart);

And my question is about the "scale" setting in the ui.Chart.image part, I don't understand what it does and therefore I don't know which value I should use.
I ran a series of comparisons varying the scale with different values (5, 25, 250, 500, 5000, 50000) but still don't understand.
So, can you advise me the "scale" value I should use?


texts: Extracting a subset of pixels from a large raster in PostGIS/ST_Clip result missing values

texts:
PostgreSQL: 12.4 (running on Amazon Aurora)
PostGIS: 3.0.2
I'm attempting to clip a much larger raster (global meteorological data on a 0.25 degree resolution grid, split into 64x64 tiles, with 9 bands for different altitudes) down to a subset around the point I'm interested in.
Beyond that, for context, I want to obtain a Binlinearly interpolated point from the raster. To do that, I've been using ST_SnapToGrid to manipulate the pixel centroids to align with my desired point and interpolate the data. My intention with the clipping is that rather than manipulating the entire raster, I'm only processing the part that's relevant, which I hope will improve performance.
Data has been ingested from a TIFF file using raster2pgsql and piped to pgsql:
raster2pgsql -s 4236 -I -a -M data_20200513_12z_t06.tiff -F -t 64x64 public.data_table -Y | psql --dbname $DB_NAME --host $DB_SERVER --username $DB_USERNAME -q

I'm running this query to get my subset of the original raster:
WITH myquery AS (
    SELECT
        ST_GeomFromText('POINT (148.1484375 -89.630)', 4236) as geom,
        'data_20200513_12z_t06.tiff' as filename
),
subset_raster AS (
    SELECT
        ST_Clip(rast, ST_Buffer(geom, ST_ScaleX(rast)), ST_BandNoDataValue(rast)) as clipped_rast, data_table.filename, rast as original_rast, geom as query_geom
    FROM
        data_table
    INNER JOIN
        myquery ON myquery.filename=data_table.filename
    WHERE
        ST_Intersects(myquery.geom, rast)
)
SELECT (ST_DumpValues(clipped_rast)).* FROM subset_raster

The clip seems to work, as it creates a new 3x3 pixel raster based on the expected coordinates SELECT (ST_MetaData(clipped_rast)).* FROM subset_raster:
147.875 -89.375 3   3   0.25    -0.25   0   0   4236    9

But there's no data in the raster SELECT (ST_DumpValues(clipped_rast)).* FROM subset_raster:
1   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
2   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
3   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
4   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
5   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
6   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
7   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
8   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
9   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}

So what am I missing? Is this the intended behaviour?
How should I query the data if I want to constrain the raster to minimise unnecessary processing?


texts: Saving Geo-Tif from QGIS creates black TIF?

texts:
I am attempting to take a Tif file that I received from a company and convert it to a geo-tif using QGIS.  The original file(s) I received were the *.tif, *.tfw, and *.prj.  I'm fairly new to working with geo-tifs, but have learned that the *.tfw is the world file that contains the points to line up the rasters and the *.prj holds the projection information.  So I was able to import the Tif into QGIS, but it was not in the projection I wanted.  After some research, I was able to convert the Tif to the right projection and I was ready to save it as a Geo-Tif.

However when I saved the Tif, the *.tif file I get was a solid black picture, not the colored map that I originally imported.  I can load the original TIF in GIMP with no problems, but when I try to load the TIF exported from QGIS, I get error messages "Unknown field with tag 33550 (0x830e).  I assume GIMP doesn't know what to do with the Geo-tag information and is erroring.  However, when it does load the tif, it is a solid black map.  It has lost all the colors from the original.

I've searched the site, but have not been able to find anything on this.  Any ideas where I could be going wrong?  The original files can be found at: Original TIF files zipped


texts: WFS filter query does not return points (selected from polygon)

texts:
I'm a newbie to OpenLayers (and JavaScript for that matter) and hoping you can offer some advice. I have been referring to examples http://openlayers.org/dev/examples/wfs-filter.html and http://osgeo-org.1560.x6.nabble.com/Select-WFS-Features-using-intersect-of-polygon-td4467419.html but have not had any luck getting any features returned from a spatial filter (a user-drawn polygon). The server response is as follows, and I'm not sure what it means:

<ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>
        org.xml.sax.SAXParseException; lineNumber: 1; columnNumber: 294;
        Element or attribute do not match QName production: QName::=(NCName':')?NCName.
        Element or attribute do not match QName production: QName::=(NCName':')?NCName.
    </ows:ExceptionText>
</ows:Exception>

Anyway, my code is below. Please let me know if there is something I am missing.

var wfsPIF = new OpenLayers.Layer.Vector("PIFs", {
    strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        url: "http://new-raptor.weatherfusion.com/geoserver/wf_data/wfs",
        featureType: "erie_pif_wfs",
        featurePrefix:"",
        geometryName: "geom"
    }),
    projection: "EPSG:4326",
    styleMap: new OpenLayers.StyleMap({
       pointRadius: 2,
       fillColor: "red"
    })
});
map.addLayer(wfsPIF);

var drawControl, polygonLayer;

polygonLayer = new OpenLayers.Layer.Vector("Query Polygon");
map.addLayer(polygonLayer);
drawControl = new OpenLayers.Control.DrawFeature(polygonLayer,OpenLayers.Handler.Polygon);
map.addControl(drawControl);
drawControl.activate();

polygonLayer.events.on({
    featureadded: function(event) {
        var geometry = event.feature.geometry;
        wfsPIF.filter = new OpenLayers.Filter.Spatial({
            type: OpenLayers.Filter.Spatial.WITHIN,
            value: event.feature.geometry
        });
        wfsPIF.refresh({force: true});
        return false;
    }
});


texts: Landsat ML classification: Apply accuracy measures to scenes without reference data?

texts:
I applied supervised maximum likelihood classification to a yearly stack of Landsat images (resolution 30m). The land cover classes are related to coniferous forest, i.e., forest, clear-cut, fire, bark beetle. I had only few years of available reference aerial photos (resolution 0.5m).
Aerial photos do not fully cover the extent of Landsat images (AOI). To assess the accuracy of the classification for single year scene, I extracted a overlapping part between the extent of Landsat and aerial photos:

I calculated the measures of accuracy for the extracted part of fully classified Landsat by set of stratified sampling points per mapped class.

I followed the methods published here:

Olofsson, P., Foody, G.M., Herold, M., Stehman, S. V., Woodcock, C.E., Wulder, M.A., 2014. Good practices for estimating area and assessing accuracy of land change. Remote Sens. Environ. 148, 42–57. doi:10.1016/j.rse.2014.02.015 and
Olofsson, P., Foody, G.M., Stehman, S. V., Woodcock, C.E., 2013. Making better use of accuracy data in land change studies: Estimating accuracy and area and quantifying uncertainty using stratified estimation. Remote Sens. Environ. 129, 122–131. doi:10.1016/j.rse.2012.10.031

to assess the accuracy of the classification.

Great R code summarizing this can be found here:
https://github.com/openforis/accuracy-assessment/blob/master/error_matrix_analysis.R

The error matrix, by this method, is not based on the simple sample counts, but on the estimated proportions of the area. Including the proportion of the mapped classes in accuracy assessment, it is possible to quantify the uncertainty attributable to variability in the sampling. Simply, to adjust the mapped area to map classification error. The error-adjusted area estimates should be accompagnied by confidence interval.

My problem:

I classified 7 Landsat images. I calculated accuracy for extracted parts of them for two years, because of no availability of another reference data sets.

My question:

Can I apply the accuracy measures from the years with reference photos to the years without reference photos? (i.e. overall accuracy is 90% in year 1, I assume the same accuracy is in year 2, 3, ...4?)


texts:
I'm working with code from the Disconnected Islands plugin in QGIS x64 and Python 3. Normally, the plugin uses networkx to create a geometric network and ends by writing the FID of features with no connections to a CSV and their corresponding network group ID.

I'm not familiar with network theory/geometric networks, so I'm not sure if this is possible: is there a way to return the fid of segments with only one connection? If so, which part of the code should I change?

# get the network
layer = iface.activeLayer()
G = nx.Graph()

# construct graph
for feat in layer.getFeatures():
    line = feat.geometry().asPolyline()
    for i in range(len(line)-1):
        G.add_edges_from([((line[i][0], line[i][1]), (line[i+1][0], line[i+1][1]),
                          {'fid': feat.id()})])

# evaluate on connected components
connected_components = list(nx.connected_component_subgraphs(G))

# gather edges and components to which they belong
fid_comp = {}
for i, graph in enumerate(connected_components):
   for edge in graph.edges_iter(data=True):
       fid_comp[edge[2].get('fid', None)] = i


texts:
Getting error from Firebug on running this code:

  TypeError: f is null  (ol.js 409)

The map is showing up properly. After adding the line : this.map.addLayer(vectorLayer); I get the mentioned error. :-(

var projection_lan = 'EPSG:4326';
var projection_lat = 'EPSG:3857';
var marker_icon = 'resources/img/green_marker.png';
var map = null;

var pos = ol.proj.transform([ 12.102727890014648, 49.017269896565764 ],projection_lan, projection_lat);

function init() {

    var localTiles = new ol.layer.Tile({
        source : new ol.source.XYZ({
            url : this.tiles_url
        })
    });

    this.map = new ol.Map({
        layers : [ new ol.layer.Tile({
            source : new ol.source.MapQuest({
                layer : 'sat'
            })
        }) ],
        target : 'map',
        view : new ol.View({
            center : pos,
            zoom : 10
        })
    });

    var iconFeature = new ol.Feature(
        {geometry : new ol.geom.Point(pos)
    });

    var iconStyle = new ol.style.Style({
        image : new ol.style.Icon({
            src : 'http://ol3js.org/en/master/examples/data/icon.png'
        })
    });

    iconFeature.setStyle(iconStyle);

    var vectorSource = new ol.source.Vector({
        features : [ iconFeature ]
    });

    var vectorLayer = new ol.layer.Vector({
        source : vectorSource
    });

    this.map.addLayer(vectorLayer);
}


texts: Unmanaged ArcObjects and C#/.NET Class Design Patterns

texts:
I am currently writing a set of batch processing jobs using ArcObjects in C#.NET. I noticed that a lot of ESRI examples pile the code into one static class with each member declared inside each method scope.

I am considering using either a factory pattern or a singleton pattern for processing features. The intent is to wrap the unmanaged interfaces in instance classes to compartmentalize the code a bit and help with re-usability. These classes implement the IDisposable interface to help make sure everything is properly cleaned up after use.

Considering what you know about COM objects, could you please let me know if the following class looks effective? Is my expectation of reusability reasonable considering how COM works? I.e., If I encapsulate a COM member, will it provide the same benefit as a manged type, or is this adding unnecessary overhead? As far as I know, the only managed type I am using is the Geoprocessor class, IFeatureClass is unmanaged.

public class FeatureProcessing : ProcessBase, IDisposable
{
    #region Fields

    private IFeatureClass _FeatureClassObj;

    private ITrackCancel2 _GPTrackCancel;

    private Geoprocessor _GeoProcessor;

    private bool _IsRunning = false;

    #endregion Fields

    #region Properties

    public IFeatureClass FeatureClassObj
    {
        get
        {
            return _FeatureClassObj;
        }
        private set
        {
            _FeatureClassObj = value;
        }
    }

    public bool IsRunning
    {
        get
        {
            return _IsRunning;
        }
    }

    #endregion Properties

    #region Constructor

    public FeatureProcessing(IFeatureClass feature)
    {
        _FeatureClassObj = feature;
    }

    #endregion Constructor

    #region Methods

    public void Simplify(int footTollerance = 5, SimplifyAlgorithmTypes algorythmType = SimplifyAlgorithmTypes.POINT_REMOVE)
    {
        //FeatureClass tmpFeatureClass = new FeatureClass(null);
        ESRI.ArcGIS.CartographyTools.SimplifyLine lineSimplify = new ESRI.ArcGIS.CartographyTools.SimplifyLine();
        lineSimplify.in_features = _FeatureClassObj;
        lineSimplify.out_feature_class = $"{_FeatureClassObj.AliasName}_Simplified";
        lineSimplify.tolerance = footTollerance;
        lineSimplify.algorithm = algorythmType.ToString();
        RunGeoprocessingOperation(lineSimplify);
    }

    private void RunGeoprocessingOperation(IGPProcess gpProc)
    {
        _GeoProcessor = new Geoprocessor();
        _GeoProcessor.OverwriteOutput = true;
        _GeoProcessor.ToolExecuting += Gp_ToolExecuting;
        _GeoProcessor.ToolExecuted += Gp_ToolExecuted;
        _GeoProcessor.ProgressChanged += Gp_ProgressChanged;
        IGeoProcessorResult2 gpResult = _GeoProcessor.Execute(gpProc, _GPTrackCancel) as
            IGeoProcessorResult2;
    }

    private void Gp_ToolExecuting(object sender, ToolExecutingEventArgs e)
    {
        if (e.GPResult == null) return;
        _IsRunning = true;
        IGeoProcessorResult2 result = e.GPResult as IGeoProcessorResult2;
        base.RaiseProcessEvent(ProcessEventTypes.INFO, $"GPTool: {result.Process.ToolName} started");
    }

    private void Gp_ToolExecuted(object sender, ToolExecutedEventArgs e)
    {
        _IsRunning = false;
    }

    private void Gp_ProgressChanged(object sender, ProgressChangedEventArgs e)
    {
        base.RaiseProcessEvent(ProcessEventTypes.INFO, e.Message);
    }

    public void CancelRunningTask()
    {
        if (_IsRunning != false && _GPTrackCancel != null)
        {
            _GPTrackCancel.Cancel();
            _IsRunning = false;
        }
    }

    public void Dispose()
    {
        CancelRunningTask();
        ReleaseCOMObj();
    }

    private void ReleaseCOMObj()
    {
        int refsLeft = 0;
        if (_FeatureClassObj != null)
        {
            do { refsLeft = Marshal.ReleaseComObject(_FeatureClassObj); }
            while (refsLeft > 0);
        }
    }
    #endregion Methods
}


texts: ArcMap 10.5.1 Export Feature Attribute to ASCII zeroes out certain columns

texts:
This is an export of about 20,000 points with 9 attributes to .csv.

Most fields come through perfectly fine, but two fields get filled with zeroes.  One has data type Double and the other is Long, and each holds numbers that are essentially integers (no decimal point).  Null values are not allowed.  Several other fields with the same characteristics are ok.

Is there anything I can change to get these fields to export properly, or is this a bug?

Before -- shapefile to export:

After -- exported .csv, sorted on column 7 to prove there are no other values:

Please note I selected only certain columns for export and did not bring along the field names.  Field 7 is PlotID and Field 9 is inner.  The duplicated columns are a little mysterious -- likely they were there but hidden.


texts:
I got a layer with over 4 million points. The distance between the points are 1 meter.

How can i reduce the data so that no points have less than e.g. 10 meters between them?


texts:
I am working with MODIS 1b Data in R, and converted the HDF files into tiffs. They load and plot correctly.

I am now trying to convert the values of corrected_counts_scale and corrected_counts_offset to dn** as described in the User Guide.

This question describes how to read this information using the CMD point. Is there a way to do this using R?

** edit - to be clear that I am not trying to use the Tiff files directly to access the attributes.


texts:
I have a point in the sea laguna and want to measure distance from the point to the closest shoreline within specified interval: every 10 degrees - new measurement.

How can I do it automatically (cause I have quite a lot points)?

I am using ArcGIS Desktop 10.5.1 with Advanced license.


texts: Select multiple choice of marker - leaflet draw angular 2

texts:
I want to create a function that can select multiple choice of marker like Mapmarker National Geographic using a leaflet draw component. I already created a toolbar that the user can draw polylines, polygons, circles, rectangles and put a default icon marker.

I have already done research into how to select multiple choice of markers but I have not found a suitable solution for my problem. Can somebody help me to create this function?

Here is my code:

import {Component, ViewEncapsulation} from '@angular/core';
import {MapService} from "./../../shared/leaflet/map.service";

import {GeocodingService} from "./../../shared/leaflet/geocoding.service";

import 'leaflet-map';

@Component({
  encapsulation: ViewEncapsulation.None,
  selector: 'maps',
  templateUrl: './app/components/maps/maps.component.html',
  styleUrls: [ './app/components/maps/maps.components.css' ],
  providers: [MapService, GeocodingService]
})

export class MapsComponent {

 constructor(private mapService: MapService, private geocoder: GeocodingService) {
  }

ngOnInit(): void {
    console.log('map component');
  }

  ngAfterViewInit() {

let map = L.map("map", {
      zoomControl: false,
      center: L.latLng(40.731253, -73.996139),
      zoom: 12,
      minZoom: 4,
      maxZoom: 19,
      layers: [this.mapService.baseMaps.OpenStreetMap]
  });

L.control.zoom({ position: "topright" }).addTo(map);
L.control.scale().addTo(map);

this.mapService.map = map;
    this.geocoder.getCurrentLocation()
        .subscribe(
            location => map.panTo([location.latitude, location.longitude]),
            err => console.error(err)
        );

var drawnItems = new L.FeatureGroup();
map.addLayer(drawnItems);
var drawControl = new L.Control.Draw({
    edit: {
       featureGroup: drawnItems,
    }
});
map.addControl(drawControl);

map.on(L.Draw.Event.CREATED, function (e) {
    var layer = e.layer;

         drawnItems.addLayer(layer);
     });
   }
}


texts: Adjust GeoServer/Geotools Heatmap code, to show it as tiled WMS

texts:
I try to change the GeoTools HeatmapProcess.java and HeatmapSurface.java classes (https://github.com/geotools/geotools/tree/master/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector), so that the resulting Heatmap can be displayed as a tiled WMS request without the cut off edges:

Anybody has an hint, what I have to change in the code to enable this?

My idea was the problem lies in the coordinate -> grid transformation to calculate the kernel density grid. I tried to increase the buffer radius of the grid calculation for each tile so that all points are included which are within the Kernel Density radius.

So far I tried to change the HeatmapSurface.java class by increasing the grid size by 100, increase the point offset by 50 and extract the grid respectively just for testing purposes but it doesn't have the intended effect:

private void init() {

    gridTrans = new GridTransform(srcEnv, xSize, ySize);

    gridTrans.setClamp(false);
    gridTrans2.setClamp(false);

    int xSizeExp = (xSize + 2 * kernelRadiusGrid) + 100;
    int ySizeExp = (ySize + 2 * kernelRadiusGrid) + 100;

    grid = new float[xSizeExp][ySizeExp];

}
public void addPoint(double x, double y, double value)
{

    int gi = gridTrans.i(x) + kernelRadiusGrid + 50;
    int gj = gridTrans.j(y) + kernelRadiusGrid + 50;

    if (gi < 0 || gi > grid.length || gj < 0 || gj > grid[0].length)
        return;

    grid[gi][gj] += value;

public float[][] computeSurface() {

    computeHeatmap(grid, kernelRadiusGrid);

    float[][] gridOut = extractGrid(grid, kernelRadiusGrid + 50, kernelRadiusGrid + 50, xSize, ySize);

    return gridOut;
}


texts: MapInfo v QGIS British National Grid Projection

texts:
I have a question about the British National Grid (ESPG:27700) projection and how it's defined in both QGIS, ArcGIS and MapInfo. The go to GIS editor in my industry has been MapInfo. It seems however that some parts of the industry have been shifting to Arc and/or QGIS and I have been encountering shapefiles produced with either of these more and more frequently.

The problem lies in the fact that the hydraulic modelling software (Infoworks) that we use employs the MapXtreme engine in order to project model and GIS data in the Geoplan. The projection used for all data is the British National Grid (BNG). When a shapefile, saved in Arc or QGIS as BNG, is opened in Infoworks, the data doesn't line up with existing MapInfo TAB files which are also BNG but saved in MapInfo. I'll give an example to illustrate.

I have received a polygon shapefile, produced in ArcGIS, that I open in QGIS along with a background map TAB we already have, produced in MapInfo. All layers were saved in BNG originally. Upon opening in QGIS however, the software doesn't recognise either as being BNG and creates USER CRS for both.

The ArcGIS shapefile is imported with the following CRS:

USER:100000: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

The MapInfo TAB file is imported with:

USER:100001: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 +units=m +no_defs

BNG is defined in QGIS as:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

While BNG in MapInfo is defined as:

"British National Grid\p27700", 8, 79, 7, -2, 49, 0.9996012717, 400000, -100000

I have tried saving the shapefile to BNG in QGIS as well as the USER:100001 CRS created when importing MapInfo TAB files, however when opened in Infoworks, it still doesn't line up. One option would be to reproject all the data (MapInfo files included) in QGIS. However, then the GIS wouldn't line up with the model in the Infoworks Geoplan, seeing as Infoworks still uses the MapXtreme engine.

As an additional note, after importing MapInfo TABs into QGIS and saving simply using the USER:100001 CRS that QGIS automatically assigns, the data won't line up any more. The only way I was able to get the ArcGIS shapefile to line up is by opening it in MapInfo and saving a copy with the projection set to BNG.

The problem with this is we are investigating moving to QGIS full time to avoid MapInfo or Arc licensing fees.

Is there a way that I can define a custom CRS that would enable me to save layers in QGIS that will line up with MapInfo data in a MapXtreme environment?