texts: Setting coordinate system in CARTO Builder map?

texts:
I uploaded a polygon layer in EPSG 25830, called 'parcels'.

I know that the geometries are reprojected and the dataset has two geometry columns:

'the_geom' in EPSG 4326
'the_geom_webmercator' in EPSG 3857.

To reproject the layer I added it to a map, which only has that layer. I used CARTO ‘Builder’. In the SQL option of the layer, I executed the following sentence, according to the documentation:

SELECT
  ST_Transform(the_geom_webmercator, 25830) as the_geom_webmercator,
  cartodb_id,
  the_geom,
  area
FROM
  parcels

The SQL sentence runs, but the polygons disappear out of the map, despite I pressed the button ‘Center map on layer’, and played with the zoom. The data of the table continues appearing.

I also changed the base layer to a ‘color’ type, and also a personalized WMS. The layer polygons do not appear. When I reset the layer data to SELECT * FROM parcels the polygons appears over the personalized WMS basemap.

I also wanted to know how to request personalized WMS base maps in EPSG 25830, or any other SRC, if it is possible?


texts: How to display coordinate-system properties of raster data in FME Data Inspector?

texts:
FME Data Inspector shows detailed information of coordinate-system properties of vector data sets (marked in green). But for raster data sets FME displays in the status bar only basic information (marked in red).

How can I display detailed coordinate-system properties of raster data in FME Data Inspector?

I did not found a way to do this in google and FME-help.


texts: Separate overlapping line shapes (or circle shapes) in raster data

texts:
I have a classifier that creates a (raster) segmentation mask that looks like this (after converting Raster to Polygon):

Now I'm wondering, whether there is a method to extract actual objects from the shape of the pixel masks - like fitting lines (respectively circles) on the orange pixels (respectively green) pixels to separate obviously different objects. A desired output would look somehow like this:

Do you have any ideas on how to approach this, with QGIS, ArcGIS or if need be programming (preferably MATLAB, python or C++)?


texts:
According to Wikipedia, ASCII-art is one of the earliest techniques used for digital graphical design. An ASCII-art technique can be traced back to 1939 and computer-art pioneers have used it since 1966. Nowadays ASCII art can be found in info files, console applications, non-graphic computer terminals and etc.

Around the Web one can see various beautiful examples of ASCII art, but not as that impressive examples of maps. Do you think that ASCII-art is just far too unsuitable for cartographic purposes? If not, can you post examples of an ASCII-art map?


texts: postgis SRTM3 different raster scale in same column

texts:
SRTM3 elevation data are stored in HGT files each containing 1201x1201 values. I want to store this information in a PostGIS raster column with a single band.

To optimize my queries, I wanted to split these files in a smaller set. For example 50x50 (and some with 50x51 and 51x50 for the edge values).

When I ran the raster2pgsql command, it has generated a table with this constraint :

enforce_same_alignment_rast" CHECK (st_samealignment(rast, '010000000...'::raster))

It seems that it will fail the check if I split with my previous scheme. But 1201 is a prime number ...

So my question is this : is it possible to store different raster scale in the same raster column ? And if yes, will it have any negative impact ?


texts: ERROR 000989 Issue with Python Script and Field Calculating distance between two long lat points

texts:
very new to scripting and I have been struggling with my script for a while. I'm not sure why it doesn't work in field calculator. I am trying to calculate a list of distance between two long lat points. The following is the code that I have implemented in the Pre-Logic Script Code for FieldName1. I am currently using ArcGIS 10.1

def Distance(Latitude, Latitude2, Longitude, Longitude2)
    Longitude, Latitude, Longitude2, Latitude2 = map(math.radians, [Longitude, Latitude, Longitude2, Latitude2])
    dlong = Longitude2-Longitude
    dlat = Latitude2 - Latitude
    a = math.sin(dlat/2)**2 + math.cos(Latitude * math.cos(Latitude2) * math.sin(dlong/2)**2
    c = 2 * math.asin(math.sqrt(a))
    Distance = 6371 * C * 1000

The following are sample points for reference:
Latitude :48.5761
Longitude : -124.31111
Latitude 2: 48.57611
Longitude 2: -124.31111

I kept getting error 000989: Python syntax error: Parsing error SyntaxError:invalid syntax(line1). It also adds that the parameter are not valid.


texts: Error generating chart: no features contain non-null values of "DispAlt" - Issue creating a chart in GEE

texts:
I'm trying to create a table in my Google Earth Engine project that populates data based on where the user selects with the cursor.

I converted this code from an example in earth engine and adjusted the variables for my data.  It works when I use data GEE has already loaded into their collection (ie Tigerline data) but it doesn't work with the data I have loaded into the code editor, which is a table from an pelican tracking sensor.

var data = table

Map.addLayer(data, null, 'birds', true);

Map.style().set({cursor: 'crosshair'});
// A list of points the user has clicked on, as [lon,lat] tuples.
var selectedPoints = [];

// Returns the list of countries the user has selected.
function getSelectedBirds() {
  return data.filterBounds(ee.Geometry.MultiPoint(selectedPoints));

}

// Makes a table of the given FeatureCollection of countries by name.
function makeResultsTable(data) {

  var newtable = ui.Chart.feature.byFeature(data, 'name');
  newtable.setChartType('Table');
  newtable.setOptions({allowHtml: true, pageSize: 5});
  newtable.style().set({stretch: 'both'});
  return newtable;
}
function updateChart() {

  var chart = makeResultsTable(getSelectedBirds());
  resultsPanel.clear().add(chart);
}

function clearResults() {
  selectedPoints = [];
  Map.layers().remove(Map.layers().get(1));
  var instructionsLabel = ui.Label('Select birds to compare altitude.');
  resultsPanel.widgets().reset([instructionsLabel]);
}

// Register a click handler for the map that adds the clicked point to the
// list and updates the map overlay and chart accordingly.
function handleMapClick(location) {
  selectedPoints.push([location.lon, location.lat]);
  //updateOverlay();
  updateChart();
}
Map.onClick(handleMapClick);

// A panel containing the two buttons .
var buttonPanel = ui.Panel(
    [ui.Button('Clear results', clearResults)],
    ui.Panel.Layout.Flow('horizontal'), {margin: '0 0 0 auto', width: '500px'});

var resultsPanel = ui.Panel({style: {position: 'bottom-right'}});
Map.add(resultsPanel);
clearResults();nter code here

The function makeResultsTable() currently spits out the error: Error generating chart: No features contain non-null values of "name".

From the research I have done it seems as though my issue is that GEE is not finding my client-side data inside the function when it's called.  However I'm new to GEE and not very experienced with JavaScript so I'm having a hard time reconciling this issue.


texts: What is this PostGIS query doing to measure the distance? How distance value is differed by the_geom and the_geom_webmercator

texts:
I dont know how the distance is differed based on WGS84 and  Web Mercator.

Sample:

SELECT *,

ST_Distance(
    the_geom::geography,
    CDB_LatLng(35.68879039,
        139.74517892)::geography
  ) / 1000 AS dist,

 ST_Distance(
    ST_Transform(
      CDB_LatLng(
        35.68879039,
        139.74517892),
      3857),
    the_geom_webmercator
    ) / 1000 AS distance

FROM tblrc_geom_jp

WHERE ST_DWithin(ST_Transform(CDB_LatLng(35.68879039, 139.74517892), 3857),the_geom_webmercator,1000 / cos(35.68879039 * pi()/180))

Results:

Record 1:
 .........
dist: 0.89083829568,
distance: 1.09980160229206

Record2:
..........
dist: 0.46739726238,
distance: 0.575471921344686

Why there is a difference in distance? My assumption is I should get the same distance value for dist and distance column.


texts: Reclassifying shapefile using ArcGIS Desktop?

texts:
So I have a shapefile consisting of 1000 polygon records. I would want to reclassify one field "Value" into a new blank field named "New Value". As an example:

record = 1; value = 25.6
record = 2; value = 35.4
record = 3; value = 14.1
record = 4; value = 44.6
record = 5; value = 8.7
record = 6; value = 2.1
record = 7; value = 10.4
record = 8; value = 18.3
record = 9; value = 30.9
record = 10; value = 28.5

I want all records with records from 0 to 10 as new field with label "rice".  Then 11 to 15 as "wheat", 15 to 25 as "corn", and the remaining unclassified records as "others".

To do this with ArcGIS Desktop, do I use the Field Calculator?


texts: Does the Bing Maps API have a distance to coast service?

texts: Does the Bing API contain a distance to coast service from an address? All we need to do is calculate the distance to a major coastline (Atlantic or Gulf of Mexico).


texts: Dijkstra’s graph returns infinite cost in QGIS?

texts:
I borrowed some code from the Network analysis library, modified as appropriate. The goal is a program that measures the distance between a location and a train station.

Eventually, this will loop through a series of addresses and stations, but I am testing it on a single pair for now.

The problem is, the cost returns infinite instead of a number.

This is the first time actually using any of this code, so I am not sure where to go from here.

The addresses come from an MMQGIS geocoded file, the stations come from a premade file.

 graphLayer =
 QgsMapLayerRegistry.instance().mapLayersByName('tl_2016_17031_addrfeat')[0]
 director = QgsLineVectorLayerDirector(graphLayer, -1, '', '', '', 3)
 properter = QgsDistanceArcProperter()
 director.addProperter(properter)
 builder = QgsGraphBuilder(graphLayer.crs())

 stations =QgsMapLayerRegistry.instance().mapLayersByName('CTA_Stations')[0]
 addresses = QgsMapLayerRegistry.instance().mapLayersByName('output.shp')[0]

 stationID = QgsFeatureRequest().setFilterFid(1)
 addressID = QgsFeatureRequest().setFilterFid(1)

 stationStation = next(stations.getFeatures(stationID))
 pStart = QgsPoint(stationStation.attribute('POINT_X'),
 stationStation.attribute('POINT_Y'))
 addressAddress = next(addresses.getFeatures(addressID))
 pStop = QgsPoint(addressAddress.attribute('Longitude'),
 addressAddress.attribute('Latitude'))

 #make list of points for all points of interest
 tiedPoints = director.makeGraph(builder, [pStart, pStop])
 graph = builder.graph()

 tStart = tiedPoints[0]
 tStop = tiedPoints[1]

 idStart = graph.findVertex(tStart)
 idStop = graph.findVertex(tStop)

 (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
 print cost[graph.findVertex(pStop)]

This is where it shows infinite.  When I get the code working, I will use the cost for further work.


texts:
I have a river water lines shapefile. I have a small clipped raster of ~7000 pixels of a place in India. I want the distance from each pixel of the raster to the nearest river line.

I am new to ArcGIS. I read something about Euclidean distance tool but am uncertain how to get the distance in an Excel sheet.


texts: Appropriate reference system for small-scale mapping of the Middle East

texts:
I want to create a thematic map of the Middle East that displays certain historic places and regions on a small-scale map (between 1:6.000.0000 and 1:12.000.000). The map has to center the boundary region between Iran, Afghanistan and Turkmenistan (about 60°E, 36°N) and possibly surrounding countries as well. It should have sufficient quality in order illustrate a scientific publication, i.e. thesis.
I am not sure about which coordinate reference system to use. UTM seems improper, since my region includes more than one UTM zone. My data is based on the WGS84 datum. So far I figured out that Lambert Conformal Conic projection may be the best, but even so I am not sure about which system to use.
Any suggestions?
EDIT (2020): I think that it was quite inappropriate to mark this question as "already answered by" a generic technical question about CRS definition in GIS, as I was trying to raise a cartographic discussion about a very specific region and use case at this time.


texts:
Simple calculation that for some reason I'm getting very wrong:

def convert_to_vector(r, lat, long):
    x = r*math.cos(lat)*math.cos(long)
    y = r*math.cos(lat)*math.sin(long)
    z = r*math.sin(lat)
    return x, y, z

def convert_to_lat_long(x, y, z):
    r = math.sqrt(math.pow(x, 2)+math.pow(y, 2)+math.pow(z, 2))
    lat = math.acos(z/r)*180/math.pi
    long = math.atan2(x, y)*180/math.pi
    return lat, long

def convert_to_dec_degrees(deg, minutes, sec, direction):
    dec = deg + minutes/60 + sec/3600
    if direction == 'S' or direction == 'W':
        dec *= -1
    return dec

kansas_city = [convert_to_dec_degrees(39, 5, 59, 'N'), convert_to_dec_degrees(94, 34, 41, 'W')]
vector1 = convert_to_vector(1, kansas_city[0], kansas_city[1])
k_city_latlon = convert_to_lat_long(vector1[0], vector1[1], vector1[2])
print(kansas_city)
print([k_city_latlon[0], k_city_latlon[1]])

I can't tell which function is giving me the trouble, but I suspect its the conversion back to lat/lon.

The results print to console are:

[39.099722222222226, -94.57805555555555]
[9.7509365327902, 108.92341788716132]


texts:
After upgrading from 3.04 to 3.10 my QGIS projects keep reporting "No transform available between EPSG:2100 and Unknown CRS" although the project settings dictate that EPSG:2100 should be used.
After spending a fair amount of time reading QGIS Github issues and suggestions such as this I cleaned my installation and moved to the latest ubuntugis version with GDAL 3.1.3 and PROJ 7.1.1 but the issue persists.

The error message appears whenever I try to use my raster to generate a new 3D view.
Also, whenever a new layer is created (eg a new Visibility Analysis) it gets a different CRS than the project default

It feels like something is not initialized correctly but it's been extremely hard to figure out what.


texts: QGIS: ESRI shapefile .shp no longer appears in list of file types under Add Vector Layer

texts:
I recently installed QGIS 2.18.4 after using ArcMap for many years. I want to add several existing ESRI shapefiles (.shp) created in Arc, to a QGIS project. To add a new layer in QGIS, I select the "Add Vector Layer" button and then click "Browse" to go find the layer. Then I navigate to the proper folder. However, the drop down menu at the lower right of the "Open an OGR Supported Vector Layer" window does not include ESRI Shapefile on the list. The first one on the list is GeoJSON. Everything after "G" appears to be okay (all the way down to "X-Plane/Flightgear") but any file formats from A through F have been cut from the list. Since "ESRI" starts with "E" it does not show up on the list. Strange thing is when I first installed QGIS, ESRI Shapefile appeared on this list. So I must have accidently done something to change how that list is populated and I want to know how to change it back.

If I type in the name of the ESRI shapefile that I know is in the folder, it and all other files starting with those letters appears in the prompt window that appears. And QGIS will correctly load the Shapefile. However, I want to be able to see all those ESRI Shapefiles without having to know and type in the first letter so it will appear in the prompt. I want it to appear in the file type drop down (or drop up?) menu.


texts: Point feature class from excel Spreadsheet lat long

texts:
I imported an excel spreadsheet of a couple hundred lat long coordinates that I wanted to become points. The majority of them were projected correctly except for a few were in the wrong hemisphere. I went in and manually added a (-) in front of the longitude to bring the points to the western hemisphere (it seems that the (-) was lost in translation or my file had errors).

The problem is that the attribute table now has the correct coordinates but the points are still in the eastern hemisphere.

I manually changed the coordinates in the attribute table during an editing session but the new coordinates have not been displayed on the map. Is there a way of changing the points coordinates in ArcMap without uploading an new edited excel sheet?


texts: Add only one row from a PostGIS layer to geoserver

texts:
I have a PostGIS table with 50 rows/entries.
Is it possible to add this 50 entries to 50 geoserver layers (1 layer = 1 row of the table)?


texts: How to find if map layer registry is empty?

texts:
In my script I want to open layer named my_layer if it is in the map layer registry. I use following code for it:

inputLayer = QgsMapLayerRegistry.instance().mapLayersByName("my_layer")[0]

If not found, I add it to the registry from file. However, if there are no layers in the registry, it returns an "out of bounds" error - mapLayersByName produces empty list with no first element (mapLayersByName("my_layer")[0]). How to avoid the error? Is catching it the only option, or can I detect empty registry before?

I'm new to Python and PyQGIS, so it's quite probable that I've missed something in my previous steps and that there's some more elegant way to load the layer without stacking instances of the same layer in the registry. If it's so, I'm interested both in the more elegant way to solve this task and in the way how to detect the empty map layer registry (I suspect I will need both later).


texts: Use PyQGIS to make grid on a vector file

texts:
I am new to Python and I am using the Python console in QGIS 3.10 to make the coordinate grid on a vector layer.
I had lots of vector files which are basically square shaped of 500 meters but all are of different location. I want to export this vector file but after making grid on top of that.
I want to do it in Python as I can loop it and do all the process at once.
My workflows goes like this:

Load vector layer
Open map canvas
Display vector layer in map canvas
Make coordinate grid on top of that
Export that as image

The sample image of what i want is like this:
Initial Vector layer:

After making Grids on that

Here is my code and in output it is nit creating the grids
`
fn = r"C:\Users\Asus\OneDrive\Documents\map secien\360m\parent grid\grid_no_A1.shp"
layer = iface.addVectorLayer(fn, 'grid_no_A1', 'ogr')

from qgis.PyQt import QtGui
layers = QgsProject.instance().mapLayersByName('grid_no_A1')
layer = layers[0]
project = QgsProject.instance()
manager = project.layoutManager()
layoutName = 'Layout1'
layouts_list = manager.printLayouts()

# remove any duplicate layouts
for layout in layouts_list:
    if layout.name() == layoutName:
        manager.removeLayout(layout)

layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName(layoutName)
manager.addLayout(layout)

# create map item in the layout
map = QgsLayoutItemMap(layout)
map.setRect(20, 20, 20, 20)

# set the map extent
ms = QgsMapSettings()
ms.setLayers([layer])

# set layers to be mapped
rect = QgsRectangle(ms.fullExtent())
rect.scale(1.0)
ms.setExtent(rect)
map.setExtent(rect)
map.setBackgroundColor(QColor(255, 255, 255, 0))
layout.addLayoutItem(map)
map.attemptMove(QgsLayoutPoint(5, 20, QgsUnitTypes.LayoutMillimeters))
map.attemptResize(QgsLayoutSize(180, 180, QgsUnitTypes.LayoutMillimeters))

grid = QgsLayoutItemMapGrid(layout)
out = QgsLayoutItemMap(layout)
grid = QgsLayoutItemMapGrid("hr", out)
grid.annotationEnabled()
grid.setIntervalX(0.1)
grid.setIntervalY(0.1)
grid.setAnnotationEnabled(True)

layout.addLayoutItem(out)
layout = manager.layoutByName(layoutName)
exporter = QgsLayoutExporter(layout)
fn = r'C:\Users\Asus\OneDrive\Documents\map secien\layout_export.jpeg'
exporter.exportToImage(fn, QgsLayoutExporter.ImageExportSettings())


texts: Categorizing attributes by keyword in new field with field calculator in QGIS 3

texts:
So, the geological classification shapefile of the Philippines has the 'DESCRIPT' field which contains the geological time scale of each feature.

"Oligocene-miocene" and "Oligocene-miocene (igneous Rocks)" and "Oligocene-miocene (sedimentary & Metamorphic Rocks)"
"Paleocene(?)-eocene (igneous Rocks)" and "Paleocene-eocene (sedimentary & Metamorphic Rocks)"
"Undifferentiated" and "Undifferentiated (igneous Rocks)"

However, some attributes have varying but same time scale.

How could the "same" time scales categorized by creating another field which contains the category?

Categorizing the features by spotting if a keyword is in the feature or specifying the exact name of the feature.


texts: Code Block for Calculate Field Excluding "No Data" Values in a Feature Class in ArcGIS Pro

texts:
I am trying to calculate a field (NAT_RE_MI) by simply adding three other fields (FOR_FISH, MINING, FARM_EMP) together.

Some of the fields have a no data value of -98 or -99.

What code would I use to exclude those from the calculation in the Field Calculator for ArcGIS Pro?


texts: I have time-series data (rasters) from which I need to calculate the direction of change (increase, no change or decrease) for every pixel from one time-step to the next. I further need an output where an increase is assigned a value of +1, no change is assigned a 0 and a decrease is assigned a value of -1 in the output raster. Not sure how I should go about this as I have thus far not been successful. Please help.


texts:
Can anyone possibly give me an example of how the Openlayers 3 'setGeometry' method can be used on an ol.Point feature please? I have a feature which has been created from data retrieved from a database and this effectively shows waypoints along a linestring route. When the route is modified I'm using the 'modifyend' event to fire a function which will update the underlying data. The last thing I need to do is find the corresponding marker within an ol.Collection and update it's coordinate property. Looking at the API, it would seem that I should be able to do this with:

// results is an array with the old coordinate and the new one
var originalCoordinate = results[0];
var modifiedCoordinate = results[1];

var sourceFeatures = session.markerLayer.getSource().getFeatures();

    sourceFeatures.forEach(displayFeature);

function displayFeature(item, index) {
    var featureGuid = item.get('Guid');
    console.info(featureGuid, index);

    //markerFC.push(marker);
    if (featureGuid == pointGuid) {
        var coords = item.getGeometry().getCoordinates();
        item.Longitude = nc4326[0];
        item.Latitude = nc4326[1];
        item.setGeometry(modifiedCoordinate);
        console.info(item, coords, index);
    }
}

item.setGeometry(modifiedCoordinate);

But it hasn't changing anything and I get an error in the console saying 'Uncaught Error: addEventListener and attachEvent are unavailable.'

I really don't understand why this is so complicated to do and the only other thing I can think of to do is destroy the feature and create a new one in it's place but that seems a little extreme when I'm only trying to update it's lat/lon properties.

Can anyone please advise how I can handle this?

So, to update this question, the item's geometry does seem to be updating when I set it here but the layer won't refresh and the item appears in it's original position. I've also tried the solution in this question but it's not working.


texts: Arcgis JS API 3.x Popup infotemplate (popup) with relationship attributes

texts:
I've buildt a halfly harcoded popup builder widget, but i don't know how to push fields from relationships into my popup div. I know how to do it in the ObjOriented way, but it's not right for me.
Example:

   var popupTemplate = new PopupTemplate({
     "title": "Beverly Hills Trees By Block",
     "fieldInfos": [{
       "fieldName": "Point_Count",
       "label": "Count of Points",
       "format": {
         "places": 0,
         "digitSeparator": true
       }
     }, {
       "fieldName": "relationships/0/Point_Count_COMMON",
       "label": "Sum of species tree count",
       "format": {
         "places": 0,
         "digitSeparator": true
       },

I show my scructure: first function creates the popup content div, the second calls it, if i want popup.
But i couldnt find the documentation, how to get a graphic's relationships, i can only get the attributes.
Any idea?
     function _getPopupTextContent(graphic) {
       let popupContent = ""
       popupContent += "
" +
       "
" +
       "
"
       console.log(graphic)
       //popupContent +=  "
"+graphic.attributes["relationships/1/FIELD"]+"
" not working //undefined

       for (var property in graphic.attributes){
         console.log((property.toString()).split(/[.]/))
         var dbName = (property.toString()).split(/[.]/)[1]
         if (dbName !== undefined){
          var dbTableName = dbName.split(/[_]/)[3];
         //console.log(maximoTableName)
         var fieldName = (property.toString()).split(/[.]/)[2];
         let value = graphic.attributes["DB.DB_ESRI_V_"+dbTableName+"_TT."+fieldName]

           if(fieldName == 'FIELDNAME1'&& value !== undefined){
             popupContent += "

Label	" +
             value
           }
           if(fieldName == 'FIELDNAME2'&& value !== undefined){
             popupContent += "
Label	" +
             value
           }
           if(fieldName == 'FIELDNAME3'){
             popupContent += "
Label	" +
             value
           }
           if(fieldName == 'FIELDNAME4'&& value !== undefined){
             popupContent += "
Label	" +
             value
           }
       }
           popupContent += "
" +
                           "
" +
                           "
" +
                           "

"
       return popupContent;
     }

     function _joinAndCreateInfoTemplate(){
       var mainLayer = _getWebmapOperationalLayerById(subStationLyrId);

       for(var title in popupSubIdDict){
         var popupName = _removeAccents(title).toLowerCase(); // subLayertitle
         var subId = popupSubIdDict[title]; //sub id

         mainLayer.getLayerObject().then(function (subLyrObj) {

           if (popupName == 'sublayer'){
             popupName = new InfoTemplate();
             popupName.setTitle(title);
             popupName.setContent(_getPopupTextContent);
           }

             subLyrObj.infoTemplates[subId] = {
               infoTemplate: new InfoTemplate(popupName)
             };


texts: shp2pgsql dumping double fields with comma separator instead of dot

texts:
I guess the source shapefile is the culprit, but is there any way to change to fix it through shp2pgsql?

Obviously Postgres throw an error with: INSERT INTO "public"."a" ("area","perimeter", ...) VALUES ('5282,98830007762','313,033597376014', ...);


texts: Does GDAL/OGR and WFS support 3DZM geometry?

texts: PostGIS has a support for geometry with elevation and measurement data. I wonder which GDAL/OGR formats support 3dzm data (if any)? Will I be able to load such data into to database using GDAL/OGR (so far I found only 3d support) and work with it via WFS-T?


texts: How to select features where column matches a regex and save them in a new shapefile?

texts:
I want to select features from a shapefile where the value of column matches a regex then save only the selected features in a new shapefile.(I used regex101.com to write the regex and check if it's captured in the test string. So the regex is good).
Here's what I did:

def filterByRegex(shapefile, output):

# Read the shapefile
shpfile = gpd.read_file(shapefile)

p = re.compile(
    ur'(^[\d]\shour\sparking\s[\d]am.[\d]pm\smon.fri\sexcept\sby\szone\s.*|^[\d]\shour\svisitor\slimit\s[\d]am.[\d]pm\smon.sat\sexcept\sby\szone\sk\spermit.*|^[\d]\shour\svisitor\sparking\s[\d]am.[\d]pm\smon.fri\sexcept\sby\szone.*|^visitor\sparking\stime\slimits\s[\d][\d]\smin\s.\s\[chart\]\s.\sexcept\sby\szone.*|^stickered\s-\svisitor\sparking\stime\slimits\s.\s\[chart\]\s.\sexcept\sby\szone\s[A-Z].*|^Sign.\s[A-Z][\d][\d][\d][\d].\sReserved\sfor\scarpools\swith\spermit\s[\d]am-[\d][\d]am\smon-fri.*|^Sign.\s[A-Z][\d][\d][\d][\d].\s[\d]\shour\sparking\s[\d]am-[\d]pm\smon-sat\sexcept\sby\szone\s.*|^Sign.\s[A-Z]444[\d].\sparking.\srectangle.\sgreen.white.\s.*)', re.MULTILINE | re.IGNORECASE | re.VERBOSE)

test_str = shpfile['Comments']
print test_str
res = re.findall(p, str(test_str))
print res #empty geodataframe?

# Save the selected features to a new shapefile
res.to_file(output)

filterByRegex('./Outputs/joinedInput.shp', './Outputs/output.shp')

When I run the script I get an empty geodataframe :

Name: Comments, dtype: object
[]
Process finished with exit code 0

I tried using shpfile.Comments.str.extract() and shpfile.Comments.str.contains(), but I get the same result.

Any suggestions?


texts: Is there an optimum RAM configuration to use for FME Desktop?

texts:
I have a couple of workstations with 32GB RAM where fme.exe processes often exceed 25GB, if its available.  And I'm looking to possibly add a 3rd workstation.  I noticed in an FME2014 beta build some performance degredation for fme.exe processes using large amounts of memory. Log file had tons of ResourceManager: Optimizing Memory Usage.  Please wait...  statements, so I think it was spending more time optimizing memory than it was using it.  But I'm not in a position where I can test a newer betas and different releases at the moment.

Do production releases of FME Desktop work best on a certain RAM configuration (or is there an optimum range).

Would FME benefit from 128GB of RAM? Or would there be issues?

I'm intentionally leaving versions out of this question because I'm interested in any experience with FME Desktop with a RAM configuration higher than mine.

FWIW, it's really cool that I get to ask a question about a software that can use 25GB.


texts: How to find identical street names with consecutively numbered addresses using Arcpy

texts:
I am trying to write a python script that will only select multiple rows with the same street name if they have consecutively numbered addresses and then place those street names in a list. For example:

FULLST     | Lowest Address Number (LOR) | Highest Address Number (HIR):
ASH ST     | 100                         | 200
ASH ST     | 201                         | 300
LOUIS RD   | 100                         | 200
LOUIS RD   | 300                         | 400

In this example I need to place ASH ST in the list but leave LOUIS RD out because LOUIS represents a break in addresses for a non-contiguous street.

This is the workflow I imagine but don't know the proper coding:

streetlist = []
rows = arcpy.SearchCursor("G/MSAG/MergedTable")
for row in rows:
    if row.getValue("FULLST") == rows.next(row.getValue("FULLST")) and row.getValue("HIR") == rows.next(row.getValue("LOR"+1)):
        streetlist.append(row.getValue("FULLST"))

This is the code that was successful in adding FULLST names to the streetlist:

streetlist = []
sc = arcpy.SearchCursor(STMS, None, None, "FULLST; LOR; HIR; MCN; OEN", None)
addToList=False
previousSt = None
previousHIR = None
previousMCN = None
previousOEN = None
for row in sc:
    if not previousSt:
        previousSt = row.getValue("FULLST")
        previousHIR = row.getValue("HIR")
    else:
        if previousSt == row.getValue("FULLST") and previousMCN == row.getValue("MCN") and previousOEN == row.getValue("OEN"):
            if (row.getValue("LOR")-1) == previousHIR or row.getValue("LOR") < previousHIR:
                addToList = True
                previousHIR = row.getValue("HIR")
            else:
               previousHIR = row.getValue("HIR")
        else:
            if addToList:
                streetlist.append(previousSt)
                addToList = False
            previousSt = row.getValue("FULLST")
            previousHIR = row.getValue("HIR")
            previousMCN = row.getValue("MCN")
            previousOEN = row.getValue("OEN")

Instead of adding FULLST to the streetlist I need to add a Unique ID number so that it will only select rows that need to be merged while ignoring streets that share the same FULLST but have gaps compared to the consecutive rows. If I change the append line to:

streetlist.append(previousUID)

It will only add the first row in the series to the streetlist. For example:
Smith: 1-110 ID:458 / Smith:  141-706 ID:459 / Smith: 707-3500 ID:460. Here only ID:458 is being added to the streetlist. If I change the append line to:

streetlist.append(row.getValue("UID"))

It will add the Unique ID for the row after the series. For the Smith example above, it will skip all of those and add Stone: 1-250 ID:461.
Any suggestion on how to select the rows in between these two?


texts: Local coordinate system to EPSG via affine transformation

texts:
I have a raster file with a local coordinate system. It is based on EPSG:31466 but with a specific rotation and offset.

gdalinfo output is this:

I need to transform this gtiff back to 31466, but how?

With FME I can use the Affiner transformer to do this.

Affiner parameter:

gdalinfo of output geotiff from FME:

But I need to know the workflow with gdal.

I'm using gdal 2.3.0.


texts: PostGIS: Querying raster data along a path

texts:
I've got a raster loaded into PostGIS that contains grib forecast data which I loaded using raster2pgsql using SRID 4326. To keep things simple I am not tiling the raster so I only have one row in my table. After the import, the resulting table contains two columns: rid (integer) and rast (raster).

The source grib contain data for two points of time, and covers 9 different altitudes. I'm storing two data points for each time and altitude so I've got a total of 36 bands.

How can I query this raster to get a rowset with useful column names containing all available numeric values available for a path over the ground (geometry linestring) over time?

Here's CREATE TABLE for the table that results from running raster2pgsql:

CREATE TABLE public.cat
(
    rid integer NOT NULL DEFAULT nextval('cat_rid_seq'::regclass),
    rast raster,
    CONSTRAINT cat_pkey PRIMARY KEY (rid)
);

And I can see I have 36 bands when I run this SQL:

SELECT ST_NumBands(rast) As num_bands
FROM public.cat;


texts: Combining multiple shapefiles in arcgisscripting?

texts:
I have got a lot of point shapefiles (r20010101.shp, r20010202.shp, r20020101.shp...). Each of them has the same amount of points at exactly the same location. In a column that is named after the Date e.g. "D20010101", where I have got unique values for every shapefile.

What I want to do now, is to Combine or Merge or Join all these shapefiles in a Python command.

The result should be ONE shapefile with all those point (but please NO duplicate ones!!).

The attribute table of this shapefile should look like this: FID, Shape, Id, [Value of Date 2001/01/01], [Value of Date 2001/02/02],...


texts: Getting coordinates of shortest distance between point and line MS SQL 2008

texts:
I have a table of point data (POINTS) in a SQL Server 2008 database table (No ShortestLineTo Function!). Each one has a SpatialPoint field of type geography (4326).

I have another table of line data (LINES). Each one has a set of points (geography, 4326 as well) and looks something like this:     LINESTRING (-3.6323834096187415 56.610269344608838, -3.6317436657085147 56.610169913782741, -3.631120602112166 56.610079250198808, -3.6304975456164366 56.609988583421483, -3.6302179950805238 56.609929350497836, -3.6298861194017187 56.609789936217744, -3.6292192708989894 56.609439268468428, -3.6290017792165123 56.60930733362, -3.6287490209975957 56.609112972308289)

I know how to get which line is nearest a point and the distance of it... but I can't seem to get the actual coordinates of where the two intersect (If I draw a straight line from the point to the nearest part of the line I am wanting the lat/long of where that is!)

I'm hoping there is a fairly easy way of doing this in SQL. I've trawled through the documentation but can't see anything which achieves this. http://msdn.microsoft.com/en-us/library/bb933917.aspx . I must say that the documentation isn't great.

Here is my current query which just gives me the distance.

----For testing
DECLARE @PointID INT
SET @PointID = 1
DECLARE @SearchCount INT
SET @SearchCount = 10
DECLARE @SearchDistance INT
SET @SearchDistance = 15000

SELECT TOP(@SearchCount)
    a.PointID,
    a.pointGeom.STDistance(b.lineGeom) as [DistanceFromTarget]
FROM
    LINES b,
    POINTS a
WHERE
    a.PointID = @PointID
    AND (a.pointGeom.STDistance(b.lineGeom)) < @SearchDistance
ORDER BY a.pointGeom.STDistance(b.lineGeom)


texts: Loading GeoJSON data to Leaflet overlay?

texts:
I am trying to get GeoJSON data from an API URL to an overlay.
Code;
// - GET GEOJSON DATA WITH AJAX TO OVERLAY REGION - //
// Get the data
var regions = $.ajax({
    url:"https://public.opendatasoft.com/api/records/1.0/search/?dataset=sverige-lan-counties-of-sweden&q=&lang=en&rows=37&start=0&facet=id&facet=lan_namn&facet=geo_point_2d",
    dataType: "json",
    success: console.log("Data successfully loaded!"),
    error: function (xhr) {
        alert(xhr.statusText)
    }
})

$.when(regions).done(function() {
    var regions = document.getElementById('lan_namn');
    L.geoJSON(regions, {
        pointToLayer: function (latlng) {
            return L.circleMarkers(latlng, {
                radius:6,
                opacity: .8
            })
            .bindPopup("
 " + records.fields.lan-namn + "
");
        }
    })
});
regions.addLayer(regions);

{
  "nhits": 21,
  "parameters": {},
  "records": [
{
  "datasetid": "sverige-lan-counties-of-sweden",
  "recordid": "2bf8cfbfbd622b8d44987da4a8fdc5f957d2ed01",
  "fields": {
  "geo_shape": {},
  "id": "08",
  "lan_namn": "Kalmar län",
  "geo_point_2d": [
    57.2009990015,
    16.1929043382
  ]
},
"geometry": {
  "type": "Point",
  "coordinates": [
    16.1929043382,
    57.2009990015
  ]
},
"record_timestamp": "2019-07-15T12:39:19.034000+00:00"
},
{
  "datasetid": "sverige-lan-counties-of-sweden",
  "recordid": "b55c4a362b22739db2cdc4896ce10d1be93feb7c",
  "fields": {
    "geo_shape": {
      "type": "MultiPolygon",
      "coordinates": [
      [
        [
          [
            15.401945130895863,
            62.00965766023548
          ]
      ............


texts: population density "hotspots" within ZCTA

texts:
With the available census tables (county, tract, block, zcta) and shapefiles, is it possible to  identify "hotspots" of population concentration within a ZCTA?

The ZCTA gazeteer file provides the population density per square meter and square mile, but what I've been asked to do is, essentially, to put a "stick-pin" wherever there are dense clusters of people within a ZCTA, using the census data.

Is that even possible with the census data that is available? I am new to this subject but it seems the census block is the smallest granule, is that right?  If the block's shape were to be spatially intersected with the ZCTA's shape, the ratio of the intersection in square meters to the block's total square meters could be used to estimate the population density of the intersection, assuming an even distribution of people throughout the block. But isn't that an unwarranted assumption?


texts: If there are many overlapping polygons at the location of GetFeatureInfo, how can I populate a template with information about all of them?

texts:
I have a layer with water basins, actually subbasins. Some of these subbasins (polygons) contain smaller subbasins. If I click on a smaller subbasin, a GetFeatureInfo request is sent to mapserver, and mapserver returns HTML produced by the following template:

!-- MapServer Template -->
<p>[item name="geographicalName"]</p>

While this works fine if there's only one feature (polygon) at the location, if there are two overlapping polygons it only shows information about one of them (typically the least useful, the largest one).

Is it possible to write a template that displays information about all features at the location? Some kind of loop perhaps?
If not, is it possible to tell mapserver which of the feature to show information for? It should be the smallest polygon (subbasin) at the location (because if the user wants information about the larger one, (s)he can click in the larger one but outside the smaller one).


texts: Running MosaicToNewRaster tool fails. C#

texts:
I am currently trying to run the MosaicToNewRaster tool via C#, on a set of TIFF files.  I am inputting all the required parameters but am still receiving:

ERROR 000582:  error occurred during execution.

Which leads me to think one of my parameters are off.  Here is my code below:
Geoprocessor GP = new Geoprocessor();
        try
        {
            int bands = 3;
            string imageryFolderPath = "Path\\ImageryFolder";
            string name = "texas2"

            ESRI.ArcGIS.DataManagementTools.MosaicToNewRaster arcToolEngine = new ESRI.ArcGIS.DataManagementTools.MosaicToNewRaster();
            string inputRasters = "";
            foreach (var pathToTiff in System.IO.Directory.GetFiles(imageryFolderPath))
            {
                if (System.IO.Path.GetExtension(pathToTiff) == ".tif")
                {
                    inputRasters += Path.GetFileName(pathToTiff) + ";";
                }
            }
            arcToolEngine.input_rasters = inputRasters;
            arcToolEngine.number_of_bands = bands;
            arcToolEngine.output_location = imageryFolderPath;
            arcToolEngine.raster_dataset_name_with_extension = name + ".jpg";
            GP.Execute(arcToolEngine, null);
        }
        catch (Exception e)
        {
            string error = "";
            for (int i = 0; i < GP.MessageCount; i++)
            {
                System.Diagnostics.Debug.WriteLine(GP.GetMessage(i));
                error += GP.GetMessage(i) + "\r\n";
            }
        }

When the tool would fail, I would try running it again with bands = 1 and got the same error.  int bands, string imageryFolderPath, and string name, are all dummy variables in this example.  Are there any other parameters that I need to add to make this work?  Or is there a rule I am missing with my current parameters?


texts:
I want GDAL translate in Pythonwin to process some GeoPDFs to GeoTiffs. I can't call up GDAL module in Pywin. I used the OSGeo4 installer for original GDAL install and can leverage GDAL properly from the OSGeo4W.bat file.

My Python and Pywin are from an install of Arc 10.2.2.  I created a PythonPath User variable and added C:\Python27\ArcGIS10.2;C:\Python27\ArcGIS10.2\Lib;C:\Python27\ArcGIS10.2\DLLs;C:\Python27\ArcGIS10.2\Lib\lib-tk;C:\OSGeo4W;C:\OSGeo4W\lib;C:\OSGeo4W\lib\tk8.5;C:\OSGeo4W\bin

After a reboot I still can't call up the OSgeo/GDAL module.

While I would like to use GDAL commands in Pythonwin, I have also tried the work around

os.system("gdal_translate -of GTiff " + sourcefile + " " + destinationfile)

which is launching cmd prompts but not processing the files.

Can someone help with sorting out the GDAL module in Pywin or to clue me in on why the os.system command is not calling up GDAL translate.


texts: Load GeoJSON and then export it as an image programatically in QGIS with Python?

texts:
I have around 6000 geoJson files that I would like to import each, change the feature colors and background color, and export as just images (png/jpeg/etc doesn't matter).

So some pseudocode

for geoJsonFile in geoJsonDir
    loadIntoQGIS(geoJsonFile)
    changeSomeFeatures()
    exportAsPNG()

I found the exportAsPNG() but I get an AttributeError that QgsApplication does not have exportItemsPNG()...so I guess it's not actually in the standard API?

Here is my basic code I have right now to try and load one geoJSON and save it...

from qgis.core import *

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

vlayer = QgsVectorLayer("my_geojson.geojson","my_geojson_vector","ogr")
qgs.exportItemsPNG()

qgs.exitQgis()

print("I guess it worked?")


texts: How do I make atlas with same pages but different layer sets?

texts:
I have a project in QGIS with multiple layers that I have to render separately. E.g. parks and distances to them (their buffers), travel time isochrons (interpolation with colors), density heatmap, etc.

I have several areas of interest, and for each of these I want to generate N pages. Each page is not a single layer, but a combo (heatmap comes with contour lines for better readability).

Right now I can render one atlas for just one setup (a page is shown below). How can I scale this to several layer setups?


texts: Cut out polygons that at least partially fall in another polygon?

texts:
I have a large polygon shapefile that describes a river catchment. I also have another polygon shapefile that contains a number of smaller polygons that represent districts.

I would like to create a new shapefile, containing only those wards that are at least partially in the large polygon (the river catchment).

I think that clipping will get rid of those parts of the districts outside the river basin. However, I want to retain those parts of the districts that are partially inside the basin.

What function can I use?


texts: Converting raster to im object for point process model covariate in R?

texts:
I am trying to convert a raster object to an .im object for use with a point process model in the spatstat package in R.
I begin by creating the raster from a tiff file using the raster() package. No problems there. I then proceed by cropping
the raster according to a given extent. Again, no problem there. I then specify a spatial window (owin) defined using the same extent.
Still no problems. When I then proceed to the final step of converting the raster to the im object using as.im(), the function runs and
the new im object is created, but it has somehow lost the pixel information that was contained in the original raster such that each
pixel now has the same value in the im object.

The date file used is here:
https://www.dropbox.com/s/n67djm3n0tfa6sx/MGVF_2001_30_arc_sec.tif?dl=0

I also have a pre-cropped file, which is much smaller in size, here:
https://www.dropbox.com/s/co2rsb4vbgn5kog/mgvf.2001.africa.tif?dl=0

And the R code is as follows:

library(raster)
library(spatstat)

# First set the geographic extent we'll be using
e <- extent(-20, 60, -40, 35)

# Then read in the Maximum Green Vegetation Fraction tiff and crop it
mgvf <- raster("MGVF_2001_30_arc_sec.tif")
mgvf.2001.africa <- crop(mgvf, e)

# Now let's create a window for in spatstat
SP.win <- as(e, "SpatialPolygons")
W <- as(SP.win, "owin")

# Finally, we create the .im object
mgvf.img <- as.im(X = "mgvf.2001.africa", W = W)

# Notice, there are no errors thrown. However, compare the plots below and see the loss of information:
plot(mgvf.2001.africa)
plot(mgvf.img)

Incidentally, I have tried the above as shown as well as trying to replace the NAs in the raster prior to converting to im. The result is the same.


texts: Does r.stream.basin still exist? GRASS-QGIS2.18.7

texts:
So I'm following steps in another topic that seems to date to January of this year and is exactly what I want to do, here's the link: Get drainage basin of a polygon  and at one of the last steps in Mr. Micha's explanation (who has helped me a lot this past week) in which he suggests to run r.stream.basins using the flow direction map as input and intersection points as parameter, I can't find any algorithm asking for such a parameter. I have everything that is needed (Flow direction as well as Intersection points) and have been browsing for a few hours to no avail. I'm wondering if it's my QGIS version or simply GRASS GIS 7 not having such a function.

Has it changed name and description?

Does it still exist?


texts: QGIS 3 - shortest path routing for multiple pairs of origins and destinations over a road network

texts:
I'm trying to map routes from an O-D survey with about 450,000 individual trips with different start and end points defined by lat, long. I want to estimate total VMT for the area by extrapolating from that. I'd been hoping to use the (really great) QNEAT3 plug-in to calculate the network cost of each route.   But I've run into a problem.  QNEAT3 generates a complete O-D matrix, measuring the network distance between each and every origin point and every other destination point.  That's not a problem at all when dealing with a smaller data set.  Since I've assigned a shared trip ID to each origin and destination pair, I'm able to just select the shortest path cost calculations I'm interested in by selecting rows where the route ID is the same for both.  Then I can just ignore the "superfluous" data.  The real problem is when the data set is so much bigger. Even parceling out 20,000 pairs at a time - something I can't feasibly do with so much data to get through - the number of extraneous calculations is so big that the operation inevitably fails after about 5 or 6 hours of processing.

I've spent a ton of time on this and I'm willing to learn any new skill that I might need to do this calculation.  I'm just not sure what that might be.  I'm very familiar with QGIS. I could probably struggle through something in R if that works.  And maybe with step by step support could try to pull something off in PostgreSQL or Python.  But I'm feeling pretty stuck.

Does anyone have any suggestions of how to go about this?


texts:
I've written a python script (ArcGIS 10.0) for importing point objects from Excel sheets to our ArcSDE database. This works fine, but I'd like to make this script a bit more custamizable by the user and therefore want to use a python list for updating my data. I've tried the following, but unfortunately the update cursor does not work together with a python list and I'm getting the error "Field item does not exist".

import os, arcpy

fc = r"C:\temp\temp.gdb\tempFC"
rowlist = [("NAME", "Testname"),("COMMENT", "Testcomment")]

print "Normal UpdateCursor"
fcCursor = arcpy.UpdateCursor(fc)
for row in fcCursor:
    row.NAME = "test 1"
    row.COMMENT = "test 1"
    fcCursor.updateRow(row)
del row
del fcCursor
print "Finished"

print "UpdateCursor with python list"
fcCursor = arcpy.UpdateCursor(fc)
for row in fcCursor:
    for item in rowlist:
        row.item[0] = item[1]
    fcCursor.updateRow(row)
del row
del fcCursor
print "Finished"

I've also tried

variable = "row.%s = %s" % (item[0], item[1])
variable

instead of

row.item[0] = item[1]

I'm not getting an error anymore, but nothing will be updated then.

Anyone another idea how to make this work?


texts: Getting surrounding pixel coordinates from ee.Geometry.Point

texts:
Not the best example, but say I have an image of Sea Surface Temperature and a point like so:

 // Load an image of sea surface temperature (SST).
var sst = ee.Image('NOAA/AVHRR_Pathfinder_V52_L3/20120802025048')
  .select('sea_surface_temperature')
  .rename('sst')
  .divide(100);

var geometry = ee.Geometry.Point([-81.2338, 42.4274])

Assuming variable "geometry" to be the center point, I would like to get the pixel coordinates 1 pixel west, northwest, north, northeast, etc. From there, I would try to get the pixel values, but I think I can figure out the rest. Essentially, I am trying to create a 3 pixel by 3 pixel window of sorts using a single center point as the input to construct said window.


texts: Execute Oracle stored procedure from ArcSDESQLExecute

texts:
I'm attempting to execute an Oracle 12c stored procedure using ArcSDESQLExecute:

import arcpy

egdb_conn = arcpy.ArcSDESQLExecute(r"Database Connections\con1.sde")

sql = '''execute update_sle (2500000, 500000)'''
egdb_return = egdb_conn.execute(sql)

print "Complete."

However, I'm getting the following error:

Message File Name   Line    Position
Traceback
    <module>    <module1>   21
    execute C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\arcobjects.py  27
AttributeError: ArcSDESQLExecute: StreamExecute ArcSDE Extended error 900 ORA-00900: invalid SQL statement

Question:

How can I execute an Oracle stored procedure by using ArcSDESQLExecute?

Notes:

The python script executes successfully on select queries. Therefore, I know the script works.
I can execute the procedure from SQL Developer successfully. Therefore, I know the procedure works.
I tried SET FEEDBACK OFF but got an error. More info here.
Related question/answer here, but I think it's for SQL Server.


texts: ST_Within and ST_Contains not working for Latitude and Longitude in PostGIS

texts:
I have a table "tbl_location" in PostGIS which has latitude and longitude, so I created a new column called "point_location" :

update "tbl_location"  set point_location =
ST_Transform(ST_SetSRID(ST_MakePoint(lng, lat), 4326), 26913);

I have another table which is available here has MultiPolygon field and I am trying to find in which polygon a location(lat/lng) belongs. For that, I wrote the following query:

SELECT a.id, b.csa2010
FROM baltimore."vs14_housing" b, baltimore."tbl_location" a
WHERE
ST_Intersects(a.point_location,  b.geom)
    or
ST_Contains(a.point_location, b.geom)
    or
ST_Within(a.point_location, b.geom)

It is not returning any records. Data type for geom column in vs14_housing is geometry (MultiPolygon). Let me know what I am doing wrong.


texts: Creating rasters with exactly 1km x 1km resolution in QGIS?

texts:
I am using QGIS 2.18.11. I want to create grids with cell size that is exactly 1km by 1km using Raster>Conversion>Rasterize(Vector to Raster) and Raster>Grid(Interpolation). My area of interest is southern Pennsylvania. I need the cell size to be exact because my variables are measured per-square-km, and I need to be able to sum over grid cells.

I have tried all sorts of different combinations of CRS for the project and for the vector shapefiles I am converting, and set the cell size at 1000 pixels by 1000 pixels, but I cant get it to give me a raster with exactly 1km by 1km resolution. When I measure the cell size for the raster that is created, it comes out at around 1.3 km, with the vertical dimension shorter than the horizontal dimension.

Curiously, I can create a vector grid that is exactly 1km by 1km, but I dont see how that helps me.

What combination of project CRS and shapefile CRS can I use?


texts: Adding a new field to table in ArcCatalog

texts:
I have a command button in ArcMap, this command button working to add a new field to a specific table in ArcCatalog.

e.g. i have a table contains many fields like ObjectID , FeatureName, ...
i want to add a new field programmatically while ArcMap is running .

Is there anyway to add a new field to ArcCatalog while the ArcMap is running.


texts:
I am using Qgis 1.8 and Python (on Ubuntu 12.10) for a couple of month now, but today i got suddenly the error

Couldn't load Python support library: Cannot load library /usr/lib/qgispython: (/usr/lib/libqgispython.so: cannot open shared object file: No such file or directory)

while starting qgis. I reinstalled Qgis but nothing happend. Is there a possibility to reinstall this library?


texts: Please clarify whether I need a Bing api key for this app

texts:
In my application I need to use Bing maps as a base map.
http://svn.openlayers.org/branches/openlayers/2.11/examples/bing.html(open in IE)

I followed this example.
For this there is no bing key is added in that example.
Only they included below .ashx file...
Can I include this .ashx (http://ecn.dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=6.2&mkt=en-us) file ?
In that .ashx file ther is credentils="some key is there" , so thats y i got the doubt.. Can I use it legally?
Please clarify my doubt ??

Any help is greatly appreciated..


texts: Is there a performance penalty when using security-enabled services in ArcGIS Server?

texts:
I enabled security on our installation of ArcGIS Server last week, and modified an existing application of ours to use the service tokens to access the services.  The application works with the new services, but the performance feels...slower.  I'm not sure how else to describe it.  Is there a performance penalty when using secure services?  Has anyone using secure services had to do any additional performance tuning after enabling security.

Thanks.


texts: How do I visualize a time and space sequence of events in ArcGIS?

texts:
I would like to visualize this scenario:

I have 5 weather stations. I would like to present a sequence of rainfall events for them. For example, rain was starts at station 2 at 6:09 and then went to station 5 after 8 minutes and ...

Please guide me on how I can visualize the aforementioned scenario as a snap shot via GIS.

Suggestions with tools description are highly appreciated.


texts:
I am using the style editor in GeoServer to make some symbols for an OpenLayers based Web GIS solution. But I face a challenge when trying to rotate some symbols. I have a Point Layer that needs to be rotated based on the values in a field named "ROTATION".

This is fairly simple and I created the following SLD that do this:

But unfortunately the values in the rotation field is not providing the correct rotation (it might be radians or something)... so I get the following result in my Web GIS (it is the small black lines):

I wish the black lines to be perpendicular to the red and blue lines. I have calculated that I will get the right rotation using the formula:
(360 - "ROTATION" + 90) / 360

At least it works when I calculate a new field on the feature class using MSSQL (where the layer is currently stored). Unfortunately I can't do the trick in MSSQL as the data are currently in production.
So my question is - how do I get this formula to work on the field ROTATION in the SLD editor.

I am using GeoServer 2.10.

I have so far read that I might need something with "ogc:Literal"Formula"/ogc:Literal", but what should I write exactly.


texts: Is it possible to add the time of editing to the table of a shapefile?

texts:
I would like to automatically track the time when a point is created or edited.
I first thought about virtual fields and the now() function. However, that field always contains the current time. So maybe it would be possible to only update the virtual field one time, but I have no clue how. Is there a way or any other possibility to automatically track the time?

In contrast to the question Making automatic updated date and time field in QGIS? or rather the answers I don't think it is reasonable to select a field or save the shapefile every time after editing a new feature. So yes the quesions may be similar, but the answers of the other question don't answer my question.


texts: Nested polygon count and shapefile query

texts:
I have attached a subset of shapefile polygons with a number of nested polygons. Wondering

How would I go about writing a query to figure out how many such nested polygons exists in one shapefile.
I have written a skeleton query, looking for some guidance on this. I want to compare the area of the three polygons and code it in the parent query such that I get the polygon with the largest area.

DECLARE @geoPoint geography;
DECLARE @geog1 geography;
DECLARE @geog2 geography;
DECLARE @geog3 geography;
SET @geoPoint = geography::STGeomFromText('POINT(long lat)', 4326);

----Parent query--------------
SELECT * from [GISTable]  where geog.STContains(@geoPoint) = 1

----poly big polygon --------------
SELECT @geog1 = geog FROM [GISTable] WHERE id = 1;
SELECT @geog1.STArea();

----poly medium polygon  --------------
SELECT @geog2 = geog FROM [GISTable] WHERE id = 2;
SELECT @geog2.STArea();

----poly small polygon --------------
SELECT @geog3 = geog FROM [GISTable] WHERE id = 3;
SELECT @geog3.STArea();


texts: Can caving systems be represented in GeoJSON/TopoJSON?

texts:
I've been looking over the GeoJSON/TopoJSON specs (as best I understand it TopoJSON is an extended GeoJSON standard which uses 'arcs' to combine features into topography rather than have each polygon as a discreet object, please set me straight if I've missed something core...)

3.1.1 states that:

  Altitude or elevation MAY be included as an optional third
     element.

Which leads me to believe that GeoJSON should be able to cope with representing Topography (Side note - took me too long to work out topo-graphy != topo-logy, urgh greek!). So assuming we can do mountain ranges with it with contours as features...

What about caving systems? Can we do GeoJSON/TopoJSON cave maps - guess these would need to be similar altitude based contour poly slices. But I can't seem to find any good examples - which usually means this is a stupid idea for a reason I can't find.

Anyone got an insight? Or are caves beyond the scope of these JSON standards and better represented in alternative formats?


texts: Kriging with expectation of 0, beyond the data supports

texts:
I am trying to Krige in ArcGIS from point samples of a chemical concentration.  Towards the center, most of the points are larger, and on the periphery, all of the points are 0. The generated raster looks good in the center and towards the edges, but beyond the outside points of zero, the raster values trend back up. I set the major range of the semivariogram to be 50, and the trending upwards is beyond this range of influence for any of the points that are greater than 0.  So it seems that Arc is setting an expectation somewhere around the average of all the points. Is there are way to set the expectation to 0 in ArcGIS. If not, what is a simple way to Krige in R or QGIS where you can set the expectation and the parameters of the semivariogram?

Here's a screenshot:

See how near the periphery, especially upper left and lower right corners, the value trends back up, when it should be zero.


texts: Adding python module to QGIS python interpreter on mac?

texts:
I'm on OS X El Capitan running QGIS version 2.18.2 which I installed via the Kyngchaos installers.

I'm trying to install a plugin, which requires, among other things, that Pysal needs to be installed.  However, despite installing it via pip it's still not working.

It appears that QGIS uses a different interpreter, which I get, but I can't find it or figure out how to install anything to it.


texts: What is the cause of 'tearing' of polygons (artifacts) using R, ggplot and geom_polygon?

texts:
Thanks to the answer given in this question I have been able to subset and draw a map of electoral divisions in part of the UK, in this case Pembrokeshire. The resulting data frame is large and contains Ordnance Survey data so it would be difficult to post here, but the data frame looks like this:

> str(bar)
'data.frame':   134609 obs. of  7 variables:
 $ long : num  214206 214203 214202 214198 214187 ...
 $ lat  : num  207320 207333 207339 207347 207357 ...
 $ order: int  1 2 3 4 5 6 7 8 9 10 ...
 $ hole : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ piece: Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ group: Factor w/ 82 levels "Amroth ED.1",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ id   : chr  "Amroth ED" "Amroth ED" "Amroth ED" "Amroth ED" ...

I fed the resulting data frame to ggplot using the following code:

ggplot(bar, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour = "black", fill = "grey50")

This generates the following image, which looks nice and clean.

Then I combined this with a data frame containing population data, which looks like this:

> str(mydf)
'data.frame':   60 obs. of  22 variables:
 $ ward.code  : chr  "00NSPH" "00NSPJ" "00NSPK" "00NSPL" ...
 $ id         : chr  "Amroth ED" "Burton ED" "Camrose ED" "Carew ED" ...
 $ la         : chr  "Pembrokeshire" "Pembrokeshire" "Pembrokeshire" "Pembrokeshire" ...
 $ total      : num  1237 1737 2458 1570 1976 ...
 $ age.0.4    : num  34 86 81 92 107 76 131 77 90 95 ...
 $ age.5.9    : num  45 93 83 80 138 82 111 85 132 75 ...
 $ age.10.14  : num  65 116 123 103 111 79 151 80 135 83 ...
 $ age.15.19  : num  69 90 161 126 117 93 150 87 139 103 ...
 $ age.20.24  : num  42 63 116 58 81 63 120 58 114 79 ...
 $ age.25.29  : num  46 63 73 60 86 56 90 51 108 67 ...
 $ age.30.34  : num  38 60 87 72 99 54 115 62 76 42 ...
 $ age.35.39  : num  53 105 104 82 110 81 91 76 121 82 ...
 $ age.40.44  : num  70 142 128 107 116 88 161 89 151 92 ...
 $ age.45.49  : num  71 138 172 122 128 109 192 116 190 104 ...
 $ age.50.54  : num  93 136 204 108 133 119 168 125 174 99 ...
 $ age.55.59  : num  126 129 235 125 149 108 179 137 175 106 ...
 $ age.60.64  : num  139 162 248 170 194 129 236 183 199 136 ...
 $ age.65.69  : num  110 110 205 95 129 143 172 128 167 130 ...
 $ age.70.74  : num  81 85 174 52 100 75 110 88 113 128 ...
 $ age.75.79  : num  78 54 130 58 74 70 72 68 119 114 ...
 $ age.80.84  : num  38 50 84 33 56 43 63 42 94 62 ...
 $ age.85.plus: num  39 55 50 27 48 42 36 55 85 84 ...

...using the following code:

foo <- merge(mydf, bar)

and plotted the result like this:

ggplot(foo, aes(x = long, y = lat, group = group)) +
   geom_polygon(colour = "black", fill = "grey50")

The problem is that the resulting plot has artifacts as shown in the image below:

So, the original data frame subset from the shapefile is fine, but the merged data file has 'issues'.

Q. What might be the cause of this kind of artifact? I understand that without the full code and data this is guesswork and I apologise in advance for this but the object is very large and there may also be redistribution issues. Any hints, pointers, suggestions as to where to start looking would be appreciated.


texts: postgis create 100 foot linestring from scratch in well known text

texts:
I am beginning to learn about postgis linear referencing functions and just generally working with lines.

I want to simply make a 100 foot line in 3424 (NJ stateplane)

first I tried creating the geometry without any SRID

SELECT st_length(ST_GeomFromEWKT('LINESTRING(0 0, 100 100)'))

I figured the 0 0 to 100 100 would be a 100 foot line?

this measures at 141.4 feet..why?

I then just simply added st_setsrid and still got the same measurement

SELECT st_setsrid(ST_GeomFromEWKT('LINESTRING(0 0, 100 100)'),3424) geom

what I am not understanding about constructing a linestring??


texts: How to create proj4 definition of local coordinatesystem

texts:
i have to generate a proj4 definition string from affine transformation parameters (A-F).

For example:

I have a local coordinatesystem based on EPSG:31466

The affine parameter are:

  A = 0.933176825,
  B = -0.358394718,
  C = 2557551.066,
  D = 0.358262137,
  E = 0.933522886,
  F = 5635190.282

Coordinates of 31466:

2575489,2435

5729759,4565

coordinates of local coordinatesystem:

50600.9577413499

81876.1639280757

but how to translate the affine parameter to a proj4 string? i need it for gdal and geoserver


texts: MapInfo to ArcGIS workspace conversion not retaining all map layers in TAB Reader

texts:
I've recently come across the TAB Reader extension for ArcGIS Desktop and am using it to convert from MapInfo to ArcGIS. I'm only testing the waters at the moment, as I've been tasked with finding a quicker solution (quicker than recreating from scratch) for transferring workspaces from MapInfo to ArcGIS.

I have a MapInfo workspace with 7 separate maps in. When I go to use the WOR to ArcGIS function in TAB Reader, only 4 of these maps appear in the map list. I've tried removing these 4 maps from the MapInfo workspace, to see if the remaining 3 appear (thought the function could perhaps be limited to 4 maps at a time). However, they didn't show up. These remaining 3 maps don't include anything too drastically different from the 4 which did appear so I'm at a loss as to why this is happening.

Are there any obvious reasons as to why the WOR to ArcGIS function in TAB Reader wouldn't show all maps in the map list?


texts: Iterate through shapefile and output an image for each feature using sf in R

texts:
I would like to create some statistical profiles for a range of features; specifically, 2016 Census profiles for the Census Divisions (CD) within British Columbia (BC), Canada. This would require a map showing each CD highlighted with respect to the other CDs (28 in total).

To do this using the sf package in R and the Statistics Canada Boundary files for 2016 Census Divisions (ZIP format) [~32MB]
I would follow these pseudo-code steps:

Load the source .shp
Filter on BC (PRUID = 59)
Select the first feature/CD
Plot the entire set of features (black line with white fill)
Plot (with add = TRUE) the selected CD with some colour (black line with coloured fill)
Move to/select next feature/CD and go to Step 4

How would I code this in R?

The idea would be to generalize this so one can cycle through any set of features within a shapefile and output individual features highlighted within the entire geography.

Related: How do I generate an image of each polygon in a shapefile?


texts: Showing attributes form when modifying features

texts:
QGIS has this neat feature to show up a form, allowing to set/change attribute values, when adding features to a layer, called attributes form. It is possible to show or hide it when adding new features:

However, I would like to show up this form also when modifying geometries of already existing features. For example if I move a point to another location.
Of course I could simply open the attribute table afterwards or use identify tool to show up the form after I moved it. But I am looking for a way to automatically show up the form when e.g. moving features, just like it does when adding new features.
Any ideas how (or if) this could be done? I couldn't find any settings in QGIS 3.4 nor 3.14, only the checkbox to "apply default value on update". This indeed updates values when moving a feature, but does not open the form so users could add extra infos.


texts: Normalized ancular choice (NaCH) in sDNA

texts: Based on previous experience with depthmap, I would like to inquire about a question: does sDNA has the equivalent value of NaCH (Normalized Angular Choice) or Total Depth?


texts: Lat long co-ordinates in geojson format for Alaska(smaller), bigger Guam and Puerto Rico?

texts:
How do I get the geojson format of Alaska, Guam and Puerto Rico in different sizes ? The sizes I have for Alaska is too large and too small for Guam and Puerto Rico. I'm very new to mapping technologies. So I'm not aware of any tools that will allow me to resize these polygons to put it on my map.It would be great if somebody could point me where I could find the rescaled version (in geo Json) of these states. Also it would be awesome if someone could point me to simple tools to that would allow me to resize the co-ordinates(in geojson format)I have for these states.

Update: the sizes of Guam and Puerto Rico in the current zoom level of my map is so tiny, it's virtually impossible to launch a pop up on clicking those states.


texts: Force QGIS to not use 'Apply default' when running processing tasks

texts:
I have Apply default value on update on for some fields in the attributes form but I don't want it to update it when I run Field Calculator functions on any of the fields even though a different field maybe updated.

Rather than having to go into the form and ticking this off in each field everytime I need to run an algorithm, is there a setting somewhere to ignore the Apply default when using Field Calculator or running models/processing algorithms?

Ex

or


texts:
I am very bad at Python and I can't even get this script to work and the worst about this is that its really small script.
Let's say I have a shapefile in C:\temp\test.shp and I want to list all the field names in that shapefile.
Let's say a featureclass with 4 fieldnames: TYPE, NAME, SPEED, AREA
and I want to list the fieldnames. I thought this worked but it did not:
import arcpy

fclass = "C:/Data/Municipal.gdb/Hospitals"
fieldnames = [f.name for f in arcpy.ListFields(fclass)]

What do I need to think about if I have a featureclass in .GDB or want to use a featureclass from MXD.
Does Python list in window the field names or can I get an output to like .txt I am interested in both ways so I can learn.


texts: How to get clicked coordinates on the line in OpenLayers?

texts:
In OpenLayers, we get the line feature on the clicked point by such kind of codes:

      map.on('singleclick', function(e) {
        var hit = false;
        map.forEachFeatureAtPixel(e.pixel, function() {
          hit = true;
        }, {
          hitTolerance: hitTolerance
        });
        if (hit) {
          style.getStroke().setColor('green');
          statusElement.innerHTML = '&nbsp;A feature got hit!';
        } else {
          style.getStroke().setColor('black');
          statusElement.innerHTML = '&nbsp;No feature got hit.';
        }
        feature.changed();
      });

But, with this, we can't get exact point coordinates on the line, because of hitTolerance setting's effect.
How to get nearest point coordinates on the line?


texts:
I have the following WFS request listed here: 'https://localhost:2222/geoserver/wfs?service=WFS&' + 'version=1.1.0&request=GetFeature&typename=CAICWeb:hospitals&' + '&CQL_FILTER=' + cqlfilter + '&' + 'outputFormat=application/json&srsname=' + prj + '&' + 'bbox=' + extent.join(',') + ',' + prj;
Or should it be like this
'http://localhost:2222/geoserver/CAICWeb/wfs?service=WFS&' + 'version=1.1.0' + '&request=GetFeature' + '&typename=CAICWeb:hospitals&' + '&CQL_FILTER=' + cqlfilter + '&' + 'outputFormat=application/json&srsname=EPSG:2231&' + ',EPSG:2231';
I've just been confused by this by the ordering of variables (i.e. prj, and request).


texts:
I created a large TMS (couple of gigs) using gdal2tiles, just to discover in the end that there was a hole (a gap) in the images, similar to the example bellow.

Since then I've obtained the missing image and I'd like to recreate the TMS. Is there a way to just patch the hole in the TMS folder structure? Or do I have to generate the entire thing again?

Update: The white area in the TMS is not a missing tile. I had a missing image when I created the tiles. So my tiles have a white box in them.


texts: How can I specify a spatial reference system which has a north pole different than the standard north pole?

texts:
I have a geographic spatial reference system (i.e., latitude and longitude) which has a north pole that is not the usual north pole: it's the magnetic north pole at approximately 70.30 degrees north, 109 degrees west.

Is it possible to write a PROJ.4 string or a WKT spatial reference string describing this coordinate system? If so, what would that string be?

At the level of precision of my data, the specific ellipsoid does not matter much. Even a spherical Earth would be fine.

(Related: How can I convert geomagnetic coordinates to geographic coordinates without doing the math myself?.)


texts: Comparing Coordinate Reference Systems in GeoTools

texts:
This seems like a pretty silly problem. I'm using a Maven build with GeoTools 14.0. All I want is to compare CRSs with flexibility so that it doesn't have to be exactly equal as long as all the parameters are equal.

There's a method listed for that here: http://docs.geotools.org/stable/userguide/library/referencing/compare.html as equals(CoordinateReferenceSystem, boolean)

But it doesn't actually exist in the GeoTools library or API docs. What gives? What's the alternative besides manually doing it, which is doable but just seeing how convoluted/intricate the CRS structure is in GeoTools this seems like a lot of work that should be unnecessary for a major spatial library.


texts: Choosing value of Moran's I to say existence of spatial correlation?

texts:
I am a GIS beginner. Usually we use Moran's I to evaluate spatial autocorrelation. Everyone knows the Moran's I coefficient ranges from -1 (uniform) to 1 (clustering); I = 0 means no correlation.

If I = 0.2, 0.3, 0.4 (suppose I is statistically significant; i.e. p-value is small, like p < 0.05), do we say there exists spatial correlation?

How large I coefficient we claim strong spatial correlation empirically?

Or we should proceed to run a spatial regression model, and then observe how significant on lag y or lag error coefficients to make conclusion?


texts: Exporting from QGIS to NEXTGIS gives errors?

texts:
Exporting a project from QGIS 2,18 into my new Mini Plan account at NEXTGIS and some of the jpeg files are not successful.
I get the error:

  {"exception": "ValidationError", "message": "Raster files without
  projection info are not supported.", "attr": "source", "serializer":
  "raster_layer"}

The projection info is surely there as the images/maps show up correctly within QGIS.

Can anyone explain this error?


texts: How to calculate line raster intersection statistics in ArcGIS Pro

texts:
Years ago, I lived by Hawth's Analysis Tools for spatial ecology. Since then I've been having to figure out how to duplicate the functions in ArcGIS 10.x and now ArcGIS Pro 2.5.0. In this case I am trying to duplicate the functions of the old http://www.spatialecology.com/htools/linerasterstats.php tool.

I have a 30-m raster representing dominance of a tree species (several of them actually, but let's simplify and assume it is just one) and a polyline feature class representing coldwater streams. I want to calculate the length of polyline segments within raster cells with a particular attribute (e.g., value > 0.5).

I thought I might be able to do something with converting the polylines to points and using Extract Multi Values to Points, but I couldn't figure out how to get lengths out of that. Any suggestions on that approach or a better one?


texts:
Given a PostGIS table of polygons, what is the best way to check if each geometry is inside the union (ST_Collect) of all other geometries from the same table?

I have tried the following:

WITH first AS (
  SELECT geom, row_number OVER () AS row_no
  FROM test
)
SELECT geom, row_no, (
  SELECT ST_Contains(ST_Collect(b.geom), a.geom)
  FROM first b
  WHERE b.row_no != a.row_no
) isinside
FROM first a

While this works, the subquery can be very slow with a large number of geometries. Is there a more efficient way to achieve the same result?


texts: How to generalize QGIS 3.0 Model so that I can install on multiple computers without setting up on every computer?

texts:
I am trying to create a model to work on multiple computers with little editing. I am working in QGIS 3.0. I am trying to get the model to first set the style for the layers based on a predefined qlr. Is there a way to generalize the style file path so that it seeks a folder called "Layer Styles."

For example in ArcMap i use %Scratchworkspace% for some generalization.


texts: Failure to project shapefile to Mass State Plane NAD83 (meters) SRS with ArcMap

texts:
I am working with a shapefile that I have not been able to project to the Mass State Plane NAD83 meters (i.e., EPSG:26986) SRS using ArcMap's "Project" tool. The "Project" tool reports that "no datum conversion" is available between the SRS of the shapefile and EPSG:26968.
Taking a look at the parameters of the shapefile's SRS using ArcMap, gives the following:

Projected Coordinate System: NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001
Projection: Lambert_Conformal_Conic
False_Easting: 200000.00000000
False_Northing: 750000.00000000
Central_Meridian: -71.50000000
Standard_Parallel_1: 41.71666667
Standard_Parallel_2: 42.68333333
Latitude_Of_Origin: 41.00000000
Linear Unit: Meter
Geographic Coordinate System: GCS_GRS_1980
Datum: D_GRS_1980
Prime Meridian: Greenwich
Angular Unit: Degree

This looks weird. I've never seen an SRS that looks like this one. It seems as though the shapefile has two (different) notions of which datum its SRS is using: NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001 (see first bold text, above), and D_GRS_1980 (second bold text, above.)
As-is this shapefile, or at least its geometry, is useless for my purposes. ArcMap cannot correctly overlay it (using on-the-fly projection) with other data I'm using, and thus none of the geoprocessing tools will work correctly with it (e.g., spatial join with another dataset.) When an attempt is made to overlay this data with other data in ArcMap (using its on-the-fly projection mechanism), the shapefile data is way off from "ground truth", usually by around 40 meters.
This suggests some kind of datum mismatch to me, but it appears to be one that ArcMap is not capable of reconciling.
How can I project this shapefile to some well-defined, well-known SRS?
EPSG:26986 would be ideal, but I will settle for any other SRS that allows the data to subsequently be projected to EPSG:26986.
I have tried to get a handle on the origin of this problematical shapefile, but haven't been able to turn up much. It dates back possibly as long as 10 years, and the person/people who created it may or may not be around/available to answer questions. My best understanding is that it is the result of importing some data in EPSG:26986 (specifically, the Massachusetts state Road Inventory) into a transportation demand modeling tool called TransCAD.


texts: Got Null values when updating attributes with PyQgis

texts:
I'm trying to create and update a field and ,for now, just change it's value to 2.

Here is my code:

layer = qgis.utils.iface.activeLayer()
index = layer.fieldNameIndex("P")
# if index == -1: # field doesn't exist
caps = layer.dataProvider().capabilities()
if caps & QgsVectorDataProvider.AddAttributes:
   res = layer.dataProvider().addAttributes( [ QgsField("P", \
     QVariant.Double) ] )
layer.updateFields()
for feature in layer.getFeatures():
# None of these work
    index = layer.fieldNameIndex("P")
    feature["P"] = 2
    feature.setAttribute(index,(2.0))
    feature.attributes()[index] = 2
#################

layer.commitChanges()

But all I get is a new column filled with NULL values.
I can't figure out where I'm making a mistake though.

EDIT:

It seems the actual value has been updated (last value) :

     feature.attributes()
 [28965.0, 209294.0, 43.0, 207459.0, NULL, u'HANOVER', u'AV', 3901.0, 4049.0, 3900.0, 4048.0, u'Hanover Av', u'LOCAL STR', NULL, NULL, u'LOCAL STREET', NULL, 0.0, 50.0, 50.0, 5436.0, 43.0, 25.0, 209294.0, NULL, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, NULL, 0.0, 0.0, 0.0, 0.0, 1, 2]

But in the shapefile, it still shows NULL.


texts: Are these Easting/Northing coordinates being converted correctly?

texts:
I was trying to confirm some GPS coordinates before my last trip to Joshua Tree and found this list which is similar to this list (Appendix IV), but when I try to convert the coordinates to the usual lat/lon, they don't really match up.

It seems that Easting/Northing coordinates can be defined in different ways depending on map projection type or offset origin?

I tried converters like https://www.engineeringtoolbox.com/utm-latitude-longitude-d_1370.html or http://www.zonums.com/online/coords/cotrans.php?module=14   (zone 11?)

For example, Mustang Ranch is listed at 3768017, 575615 which converts to 34.05013, -116.18072:

but is actually closer to 34.05194, -116.18145:

Headstone Rock - South Face is listed at 3760548, 578384 which converts to 33.98257, -116.15139:

but the actual location is more like 33.98439, -116.15240:

Do these seem like conversion errors or measurement errors?


texts: Extract maximum raster value within a set radius to a point

texts: I have lines (rivers) that don't exactly match a raster grid (flow accumulation). If they were exact, I could just extract the raster values to points I put along those lines. However the raster values I want are often a few pixels away. Luckily the raster value I want are much larger than the surrounding pixel values. There used to be a tool that you could use to extract the maximum pixel value within a set distance from your point... but I can't remember the name of the tool.


texts:
I'm having a problem with CRS and layers in QGIS & on Geoserver
What I'm trying to accomplish to change CRS of layers in importing.

I'm importing 'test.geojson' file in QGIS, and when I open points it says that there is WGS84 CRS set, & I'm trying it to use projected CRS since I need meters in EPSG:4326

But as soon as I change CRS with right-click SetCRS/SetLayerCRS and change it to WGS 84/ UTM zone 31N all the points disappear.

Does someone know what could be the problem?


texts: How to set the extent of a rotated data frame from a rotated rectangle using ArcObjects?

texts:
I use ArcObjects to produce a map document.

I have a rotated rectangle that I want to be the exact extent of my data frame.

So, I want to rotate the data frame with the same orientation than the rectangle, then fill the map space with the rectangle.

What I want:

I tried this code:

IPageLayout pageLayout = (IPageLayout)mapDocument.getActiveView();
IMap dataFrame = mapDocument.getMap(0);

IPolygon polygon = ...

// data frame rotation
double dataFrameRotation = -44.46035122650748; // (the polygon orientation in negative)
IActiveView activeView = (IActiveView)dataFrame;
activeView.getScreenDisplay().getDisplayTransformation().setRotation(dataFrameRotation);

// setting the extent
IEnvelope zoomToExtent = polygon.getEnvelope();
activeView.setExtent(zoomToExtent);

The problem I get is that zoomToExtent means the polygon extent in the map coordinate system (xMin, xMax, yMin, yMax), so the data frame is not enough zoomed.

Any suggestion?


texts: Recommended Coordinate Reference System (CRS) for Germany

texts:
I have WGS84 latitude and longitude for data points near Munich and I'd like to use a grid-type coordinate system, like the Ordnance Survey grids in Ireland or the UK.

What coordinate system do you recommend? Does Germany have a similar national grid?

Then a bigger question, for any location how does one find out what the recommended local grid is?


texts: Why do I get "Error: Failure during raster IO" when extracting values from my raster stack?

texts:
In R I create a raster stack using the following code

#load rasters in a file list and store it to a variable
env_data<- list.files(path="C:/R_Data/WORK/env_stack", pattern='tif$', full.names = TRUE)
#create layer stack of predictors
env_predictors<- stack(env_data, quick=TRUE)
env_predictors

This works and I don't get an error. But when I try to extract values to my point data with

pc_values<- extract(env_predictors, pc_select)

I get the following

  Error: Failure during raster IO

I created/prepared my rasters in QGIS and ArcGIS (slope, aspect, climate data,...). In ArcGIS I masked all rasters by another raster (always the same) and set this raster as snap raster in the environmental settings to get the same extent. All rasters have the same resolution of 50 meter. When I visually check this in QGIS or ArGIS it seems to be ok and all rasters are congruent.

I also excluded some of the rasters and sometimes I don't get the error.

Why do I get this error? I'm just starting to work with R

EDIT:
I could limit the error to some files (the DEM and derivatives) and recreated them. That doesn't help either.
I also got the error when I try to plot one of these rasters plot() and when I write the raster stack on disk writeRaster()


texts: Calculating mean slope of lines based on DEM?

texts:
I have a DEM and a line vector shapefile (also available in PostGIS) and I need to calculate the mean slope of every line. Direction of the slope (positive or negative) is not important.

That means I should assign the DEM raster values of the start and end vertex to two new attributes (e.g. alti_a and alti_b). Then it would be possible to use field calculator and calculate the slope.

How can I do this? Are there alternative solutions?


texts: Creating Spatial Map Series with Overview Map showing changing location of Main Map Frame in ArcGIS Pro

texts:
I am trying to create a layout with two map frames. One map frame showing details in M 1:5000 the other one showing an overview map in M 1:50.000. I would like to link both views and show the extent of the detailed map frame automatically in the overview map frame. It is a map series of more than 6000 maps.
In ArcMap it was possible to fix the scale of the overview map, and dynamically highlight the extent of the other data frame. Everything can be exported to PDF with Data Driven Pages.
How do I do this in ArcGIS Pro?


texts: Extract ASCII cells geographical position

texts:
I have a task where I need to locate the geographical coordinates from a series of cells (indexed as row/column numbers) from an ASCII raster.

How can I do that in ArcGIS?

I know the header of the raster, with the starting x, y and the cell size. Is there a formula to find obtain the coordinates from that?


texts: tool to output xy from an input xy, distant, and angle

texts:
I want to build a tool where you provide coordinates, distance, and angle for input and it outputs coordinates, no feature necessary, using the distance and angle from the input coordinates.

How do I build this?

I am using ArcMap for desktop 10.4.1


texts:
I have a working WMS service under http://localhost:8080/geoserver/nis/wms?service=WMS&version=1.1.0&request=GetMap&layers=nis%3Atest

Also I have a Node JS application with OpenLayers. In this application I have a button. It gets visible after a certain condition is met. With this button I want the user to be able to load the WMS.

This is my button:

          </button>
      <label
        className={this.state.addMessageVisible ? "sc-fakeLink" : "sc-hidden"}
        style={{ gridColumnStart: "1", gridColumnEnd: "3", gridRowStart: "4", textAlign: "-webkit-center", alignSelf: "center" }}
        onMouseUp={this.onAddWMSToMyMaps}
      >
        Show WMS
      </label>

This is triggered by the button:

  onAddWMSToMyMaps = () => {
// ADD MYMAPS
window.emitter.emit("addWMS");};

which in turn triggers this on the map canvas:

    window.emitter.addListener("addWMS",() => this.addLayerviaWMS(true));}

and after all shall trigger this one:

  addLayerviaWMS = () => {
var newlayer = 'http://localhost:8080/geoserver/nis/wms?service=WMS&version=1.1.0&request=GetMap&layers=nis%3Atest4733'
window.map.addLayer(newlayer)};

The error is the following:

TypeError: can't assign to property "ol_uid" on "http://localhost:8080/geoserver/nis/wms?service=WMS&version=1.1.0&request=GetMap&layers=nis%3Atest4733": not an object


texts:
Anyone know of another source (preferably free) for these 50k Mexico Topo Maps other than  www.omnimap and www.eastview?

There is a free WMS that has them but only for UTM 15 and 16 here:
http://antroposig.ciesas.edu.mx:8399/arcgis/services/MapotecaDigitalCIESAS01/MapServer/WMSServer

The areas we need are in UTM 12,13,and 14.

Carta topográfica 1: 50 000


texts:
I am relatively new to GIS and am using ArcGIS Pro. I am making a map of a property that includes a watercourse and an environmental setback from the watercourse. I need to add the development plans to the map, which were given to me as a PDF (not georeferenced). The CAD file is not available.
Is it possible to add the PDF to the map so I can trace or draw the polygons?
I suppose it would have to be transparent so I could line it up with the property boundaries.
Or perhaps I could geo-reference it?


texts: Leaflet make legend with different circle sizes

texts:
I would like to make the legend with different circle sizes in Leaflet.

I found one examples:

https://stackoverflow.com/questions/37446283/creating-legend-with-circles-leaflet-r

https://stackoverflow.com/questions/58505589/circles-in-legend-for-leaflet-map-with-addcirclemarkers-in-r-without-shiny

https://gist.github.com/BobTheScientist/2ea0e0258d98380337524a633039385e

But it looks like the code comes from R language instead of JavaScript, hance I have no clue how to input it into my map.

I found also a different example:

http://jsfiddle.net/nathansnider/o563bg44/5/

which is exactly what I need

Unfortunately the incorporation into my map wasn't successful. I have made something like this:

  var legend = L.control({position: 'bottomright'});
    legend.onAdd = function (map) {

     var div = L.DomUtil.create('div', 'info legend');
     grades = [0, 10, 20, 50, 100, 200, 500, 1000],
     labels = ['Amount of units'],
     from, to;
//iterate through grades and create a scaled circle and label for each
    for (var i = 0; i < grades.length; i++) {
      from = grades[i];
       to = grades[i + 1];
        labels.push(
        '  ' + from);
    }
    div.innerHTML = labels.join('
');
    return div;
 };
  legend.addTo(map);

and there is nothing on my map.

The console says, that:

Uncaught ReferenceError: from is not defined
    at NewClass.legend.onAdd ((index):709)
    at NewClass.addTo (Control.js:70)
    at (index):720

For this purpose I was searching through the web how to fix this and I found an example here:

https://tommcfarlin.com/javascript-reference-error-is-not-defined/

where as an additional part of the code I added:

 function acmeReferenceError( valueExists ) {

'use strict';

 if ( undefined === valueExists ) {

/* If this is hit, then the function is being invoked
 * before the variable in question has been defined.
 */

} else {

/* If this particular conditional is hit, then the function
 * is being invoked after the value in question has
 * been defined.
 */

 }

}

but without the result either.

COnsole states the same.

Can anyone help me?


texts:
I created a few heatmaps based on 2 different shapefiles( point layers). Then I joined the data from the 2 shp files, the joined file looks perfectly fine. I created a view heatmaps based on the new shp file, unfortunately the heatmap layers appear in the layer overview (left side of the screen) without being displayed.

Does anybody know what might be the problem?


texts: Cannot see or classify weighed overlay or weighted sum rasters using ArcMap

texts: I ran weighed overlay and weighed sum on a suitability analysis in a model. neighther one is visible when added to display. The weighed sum values are also suspect (first pic). I go to properties, symbology and when I go to classify I get Error 1 (2nd pic) encountered an improper argument -failed to classify too many unique values, then when I try to classify I get error 2 (3rd pic). I have added statistics to both, it has not helped.


texts: Openlayers 5: zoom view to features in a layer

texts:
I'm adding a number of features (linestrings and points) to a layer and I want to zoom/pan the map view to the extent of those features.

I thought I'd have to do the following:

var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM(),
            }),
    ],
    target: 'map',
    view: new ol.View({
        center: ol.proj.fromLonLat([default_center_lon, default_center_lat]),
        maxZoom: 18,
        zoom: 7
    })
});

var features = [];

// Add one or more features

var marker = new ol.Feature({
    geometry: new ol.geom.Point(ol.proj.fromLonLat([coordinatesLong, coordinatesLat])),
        });

features.push(marker);

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

map.addLayer(layer);

var layerExtent = layer.getExtent();

if (layerExtent) {
    map.getView().fit(layerExtent);
}

But the map remains at the default center and zoom level, instead of zooming/panning to the extent of the features.


texts:
I'm running geoserver in Tomcat and went through the process of installing the plugin for GDAL formats.  I specifically need NITF.

Is there a URL I can post to my geoserver that confirms that I properly installed the GDAL plugin?


texts: How can I calculate the radius or curvature of a meandering river?

texts:
I want to calculate the radius of curvature of a meandering river with ArcGIS 10.1. I have used the Planform Statistic Tool, developed by the National Center for Earth-surface Dynamics, to calculate width and radius of curvature (Rc), but according to the authors of the tool, the values of Rc show a significant local scatter in a first-order curvature estimate, and they recommended that higher order schemes be used or curvature be smoothed prior to additional analysis.
These are the results that I got:

Doing a moving average of the values will be enough for smoothing the curvature? If so, doing a 3-values moving average, for example, will be enough? these are the results:

How can I interpret those data? The program has calculated 991 values (for the whole line of the river).

I'm calculating the Rc because I want to compare the radius of curvature with the width (Rc/w) to see the efficiency of the river, and with all these values I don't know how to do that.

So, there is another way to calculate the radius of curvature of the river?
I have read the other post but I think that I can't use for this case.


texts: It is possible to create a raster map using a vector dataset in GeoServer?

texts: I need to create a map of pressure shades using a point data set having each point its value of pressure. So I need to give a color to each point based on its pressure value. But after that I need the spaces between points to be fill with the interpolation color between every four points.


texts: Anyone built a custom GPFunction in ArcGIS10 and having closing issues?

texts:
I have a custom GP function in ArcGIS 10 written in C#.  When the tool is run it takes ArcGIS about 5 minutes to close (as in the application closes but in Task manager ArcMap.exe and ArcSOCP.exe are still running).  When closed ArcMap.exe will increase memory too about 170 - 250 MB.

To ensure it wasn't my code I created a blank GP Function with no code in the Execute function and return and empty paramarray value.  Ran the tool and ArcMap still took ages too close.

I then opened ArcMap again and launched the tool but this time pressed cancel instead of go.  ArcMap closed fine.

So as soon the tool enters the execute method it seems to then take ages to close when the Execute method finishes.

If I run one of the out of the box GP tools ArcMap closes fine.

Anyone encountered the same problem?