texts: DBMS error when importing XLSX into ArcGIS Pro

texts:
I am working with an XLSX, and consistently receiving either a "underlying DBMS error", or  associated with the file.

I have worked through all the suggestions and limitations from ESRI to the best of my understanding. My field names start with a letter, have no spaces, are less than 30 characters, and only letters and underscores.  I have also made sure all  my columns were consistently formatted: date, text, or number.  Further, I have resaved using different filenames, and have copied the data to new workbooks in case of corruption.  I have also restarted PRO as well as my computer.

I have 1391 rows, and 23 columns.

Is there a limit to the number of columns an imported worksheet might be allowed or potentially another issue I have not considered?

I have noticed that every time I copy the data to a new workbook, I get the error.  When I take the new worksheet and make sure all the column formatting is consistent, I don't get the error, but the worksheet does not show up (image above).  When I take that same newly created and column-formatted workbook, copy and paste it into yet another new workbook, the error returns, until I re-format the columns, etc, etc.

Due to data sensitivity, I can't provide an image example, but if necessary I will rename things to post.

I know https://gis.stackexchange.com/questions/312470 has previously been asked here, but no response was given, and that question is now deleted.

Note: I recently imported the same spreadsheet into MAP and it works just fine.  While I'd like to exclusively use MAP as it has been so much easier to use for everything I do, I am trying to adopt ArcGIS Pro due to the phasing out of MAP support.


texts:
I am using arcGIS 10.0  SP5

I have the following code :

import arcpy

arcpy.ImportToolbox(r"C:\Documents and Settings\Derek\FRAME\FORTEST.tbx", "Derek")

arcpy.TEST_Derek()

The model is names TEST
The toolbox alias is Derek

When I try and run this I get an error

AttributeError: 'module' object has no attribute 'TEST_Derek"

What I want to do is to run the model TEST from the toolbox FORTEST


texts:
I have downloaded QGIS 1.8 and installed it on my windows 7 system. When I click on the desktop icon, nothing happens.

I am currently running QGIS 1.7.4 with no problems and the installs of previous versions have always worked with no problems. I'm not sure what's going on here. Has anyone else experienced this and have a solution?


texts: Invalid return value: published arcgis server 10.0 toolbox

texts:
I have a published toolbox that is using conditional execution in one of its models:

This works fine in ArcGIS Desktop 10.0, and each sub model\script works fine when run separately from this model via the geoprocessing service.  It appears that ArcGIS Server does not like having a nonassigned model parameter. When this model is run via a geoprocessing service, ArcGIS Server sends a Invalid return value error regarding the first dataset of the branch that does not execute.

I could move the conditional branch to the client or to a middle wrapping web service instead.  However, I would like to keep the client side execution of this as simple as possible.  Has anyone here run into this problem or know of a workaround that doesn't require a separate web service or conditional logic in the client?


texts:
Here is the error notice I receive after starting GRASS GIS 7.6. I am not sure what is the source of these errors. I had reinstalled the application but it didn't fixed these errors.

What are possible remedial measures?

and here are gisrc file contents after installation.


texts: Extracting raster classe(s) and creating new raster that is SMALLER in size using QGIS?

texts:
Extracting raster values and creating new raster in QGIS? works great but there is no reduction in the new raster file size after raster calculator is run displaying just one or two classes (values) out of 8 classes... aerial tiffs segmented (using Orfeo) to capture tree canopies in this case.

Is there a tool that will complete this same process but ALSO result in a smaller file (tiff) size?

Dealing with 100mb tiffs that need to be reduced in size since they need to be converted to a shp.

100mb is apparently way too big for QGIS's polygonize tool as it hangs up in a seemingly infinite loop.


texts: Using “Show” property in Cesium to turn on and off visibility has a strange delay

texts:
I am loading some geometries through a datasource load promise like above:

Cesium.when(Cesium.CzmlDataSource.load(environment.apiBaseURL + `/geometry/all`), result => {
    this.dataSources = result;
    this.viewer.dataSources.add(this.dataSources);
});

and I am trying to switch their visibility with this callback above:

this.visibilitySubscription = this.store.visibility.subscribe(visibility=>{
        this.dataSources.show = visibility;
    });

but unfortunately the show/hide of the objects has a weird delay and I must scroll the mouse in order to see them hide or show


texts:
What I’m looking for is a concave_hull-by-attribute type solution. Or something that results in the same effect. I have a multipart point layer that I want to group by attribute with minimum overlap (fig1). The Concave hull plugin does not work on multipart layers. It wants to group ALL the points on a layer together (fig2).  Convex plug does accept hull-by-attribute but does not closely follow the points to the degree I need (fig3). I have also tried breaking the multipart layer into separate layers using the Spit vector layer plugin but that process combines all the points to one in the attribute table. So, then can’t be used with concave-hull because there are no longer any points.

How can I make this work?


texts: Wrong pixel dimension after SRS assignment with pyqgis

texts:
I use pyQGIS (Python 2. with QGIS 2.18) to assign the same reference system to some raster files in .asc format (that came out without any SRS).
Each .asc file has 1m pixel dimension.

I use this pyqgis script to convert the raster in GeoTIFF and assigning the EPSG:3003.

It works but output files have a final resolution of 16.2 x 15.6.
How to preserve the original pixel dimension?

import os
import processing
from qgis.core import *

raster_filepath = "/home/giacomo/Desktop/Livorno/QGIS/Raster/DEM_Lidar/new/prova"
outputDir = "/home/giacomo/Desktop/Livorno/QGIS/Raster/DEM_Lidar/new/"

for i in os.listdir(raster_filepath):

 if i.endswith (".asc"):

  #name fileOut
  ii = i[:-4]
  fileOut = outputDir + "/" + ii + "_rep.tif"

  layercount = os.path.join(raster_filepath, i)
  raster_layer = QgsRasterLayer(layercount, 'raster')

  # assessing the raster dimension
  raster_extent = raster_layer.extent()
  xmin = raster_extent.xMinimum()
  xmax = raster_extent.xMaximum()
  ymin = raster_extent.yMinimum()
  ymax = raster_extent.yMaximum()

  # processing
  processing.runalg("gdalogr:translate", \
  {"INPUT":raster_layer, \
  "OUTSIZE":100, \
  "OUTSIZE_PERC":False, \
  "EXPAND":0, \
  "PROJWIN": "%f,%f,%f,%f"%(xmin,xmax,ymin,ymax), \
  "SRS":"EPSG:3003", \
  "OUTPUT":fileOut})
  print(fileOut)


texts:
I'm quite new with ArcGis 10 so I’m having difficulties with something that may be simple to some. Anyway I have a list of bores that I’ve plotted using coordinates. Now I need to draw a single contour line joining the bore points by their conductivity levels. As fair as I’ve tried researching I need to convert the shapefile of the bore details to a raster file and then create a contour line with the Spatial Analyst. Is this the correct procedure? I'm having real trouble trying to figure it all out.
I've tried creating a raster file but not sure which cell assignment type etc to choose. Every time I try creating a point raster file based on the conductivity, I don’t get all the information or there are no visible points.
Could anyone shine a light on this for me and provide some help to the steps I may need to follow. I hope I’ve supplied enough information.
Thanks


texts:
Is it possible to ask R {sf} package to use a GDAL installation I prefer?

So far I found how to check external software it uses:

library(sf)
sf_extSoftVersion()

#          GEOS           GDAL         proj.4 GDAL_with_GEOS
#       "3.6.1"        "2.2.3"        "4.9.3"         "true"

I am hoping I can choose a path to my existing GDAL installation like {gdalUtils} package does:

library(gdalUtils)
# To change the path to GDAL
gdal_setInstallation(search_path = "C:\\Program Files\\QGIS 3.4\\bin", rescan = TRUE,
                     ignore.full_scan = TRUE, verbose = FALSE)

# To check the result
getOption("gdalUtils_gdalPath")[[1]]$path
#       [1] "C:/PROGRA~1/QGIS3~1.4/bin/"
getOption("gdalUtils_gdalPath")[[1]]$version
#       version
#       "2.4.0"

If this is possible, {sf} (and {stars}) package can use GDAL 2.4.0 under the hood of QGIS 3.4.

[Background]

In my previous question R stars package and set dimensions to a NetCDF file I leaned my issue was not reproducible in other environment. Now trying to see if change of environment can improve the situation. Upgrading GDAL version is the first step.


texts: Set a custom maps rotations for specific pages in the print composer in QGIS2.14.3

texts:
I have prepared an atlas in the print composer and want to set a different rotation on some maps.

How to do that?


texts: Get random point from django PolygonField

texts:
TL,DR; I want to get a random point from a polygon (potentially) using ST_GeneratePoints.

Background

I'm making a GeoDjango webservice and have a collection of UK Postcodes with respective boundaries like so:

from django.db import models as dj_models
from django.contrib.gis.db import models as gis_models

class Postcode(gis_models.Model):
      pretty_postcode = dj_models.CharField( max_length=8 )
      coords = gis_models.PolygonField( default='POLYGON EMPTY' )

I've found a delightful little PostGis function ST_GeneratePoints, it can find me random points in my coords region.

Your Challenge

Your mission if you so choose it, is to show how to use this function from within my python django app (or suggest a better way). Ideally ending up with a function like so:

from django.contrib.gis import geos
# ... other imports ...

class Postcode(gis_models.Model):
     # ... fields ...

     def get_random_point(self):
         rand_point = # code that executes ST_GeneratePoints
                      # and returns a geos.Point instance
         return rand_point


texts:
I am trying to stack a 3 band image (2 NDVI bands and 1 MDPI band). There is some error when I am trying to visualise the MDPI band (which is appearing as masked in the Inspector) and hence the stacked image is also not appearing. Where am I going wrong?
var sent2_1=ee.ImageCollection("COPERNICUS/S2")
           .filterBounds(roi)
           .filterDate('2019-11-15','2019-11-17')
           .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 40))
           .mean();

var sent2_2=ee.ImageCollection("COPERNICUS/S2")
           .filterBounds(roi)
           .filterDate('2019-01-10','2019-01-19')
           .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 40))
           .mean();

Map.addLayer(roi);
print(sent2_2,'sent2');
Map.addLayer(sent2_2,{min:0,max:3000,bands:['B8','B4','B3']},'sent2');
//Map.centerObject(roi,13);

var sent1=ee.ImageCollection('COPERNICUS/S1_GRD')
          .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
          .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
          .filter(ee.Filter.eq('instrumentMode', 'IW'))
          .filterDate('2020-03-09','2020-03-10')
          .filterBounds(roi)
          .mean();

print(sent1,'sent1');
Map.addLayer(sent1,{min: [-25, -20, -25], max: [0, 10, 0]},'SAR');

var ndvi=function(image){
  var ndvi_2_2=image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi_2_2);
};

var ndvi_1 = ndvi(sent2_1);
print(ndvi_1,'ndvi_1');

var ndvi_2 = ndvi(sent2_2);
print(ndvi_2,'ndvi_2');

var im2 = ndvi_1.select('NDVI');
var im1 = ndvi_2.select('NDVI');

var mpdi_fun=function(image){
  var mpdi=image.normalizedDifference(['VV', 'VH']).rename('MPDI');
  return image.addBands(mpdi);
};

var mpdi = mpdi_fun(sent1);
print(mpdi,'mpdi');

var mpdi_1=mpdi.select('MDPI')
print(mpdi_1,'mpdi_only')
Map.addLayer(mpdi_1,{min:-1, max:1},'MPDI')

var stack=im1.addBands(im2).addBands(mpdi.select('MPDI'));
Map.addLayer(stack,{max:1,min:-1},'stack');
print(stack,'stack');


texts: Draw gridded lines inside a feature and generate random points

texts:
I have several vectors containing features that originated from a 30 x 30 m raster. I want to generate a grid of lines (fishnet) in all of the features that are more than 30 x 30m. As shown in the image below:

After that, I want to generate a single point per feature that would fall inside the center of the grid.


texts: Intensity Range Correction - How to choose a power value?

texts:
I'm using lasrangecorrection from the lidR package to correct my Intensity values for effects of distance from sensor to target. I'm wondering how to best choose the f parameter (exponent for ratio of range to standard range). A power of 2.3 is recommended for vegetation, however this seems to flip the trend - instead of Intensity decaying as range increases, Intensity now increases as range increases (see images below). Currently, I'm using a trial and error framework to choose the power that will flatten the relationship. Is there a more efficient way to do this?

Relationship between Uncorrected Intensity & Range:

Relationship between Corrected Intensity & Range (power = 2.3):

Relationship between Corrected Intensity & Range (power = 1.6):

EDIT
After the comment from @JRR, and examining an uncorrected intensity raster, it seems correction is not necessary.

Uncorrected Intensity Map With Estimated Sensor Positions Overlayed (mean intensity at 3m resolution):

Digital Elevation Model for Area of Uncorrected Intensity Image:

Canopy Height Model for Area of Uncorrected Intensity Image:


texts: Visualizing road Intersections & connections in ArcGIS

texts:
I am working with polylines (roads) in ArcGIS and I was wondering if there is a possibility to somehow show intersections as well as connecting parts of all the different road-segments.

Also, would it be possible to present these results visually (e.g. Points on each connection/intersection) as well as in table form (“which road has how many connections/intersections” or "which road intersects/connects with which road)?

I found a few similar questions, however the answers were not providing me with usable information.
As I am quite new to ArcGIS and am starting to get frustrated finding a way to do this.

Could the Network Analyst tool be the solution for this?


texts:
I am doing the sample example from OpenGeo WorkShop. It's about creating simple popup. I get popup but with blank no feature information? I am using geoserver with geoext.

controls.push(new OpenLayers.Control.WMSGetFeatureInfo({
                autoActivate: true,
                infoFormat: "application/vnd.ogc.gml",
                maxFeatures: 3,
                eventListeners: {
                    "getfeatureinfo": function(e) {
                        alert(e.xy);
                        var items = [];
                        Ext.each(e.features, function(feature) {
                            items.push({
                                xtype: "propertygrid",
                                title: feature.fid,
                                source: feature.attributes
                            });
                        });
                        new GeoExt.Popup({
                            title: "Info",
                            width: 200,
                            height: 200,
                            layout: "accordion",
                            map: app.mapPanel,
                            location: e.xy,
                            maximizable: true,
                            collapsible: true,
                            items: items
                        }).show();
                    }
                }
            }));


texts: Should I use NAD83 datum for surveys in North America if features are traced over Google Earth (WGS84)

texts: I know NAD83 is more accurate in North America, but the features we inspect are just traced over google earth imagery which from my understanding is referenced using WGS84. Would it be better to keep the data in the original datum or to transform it into NAD83? Would any problems arise if we projected the data to NAD83 and then chose the correct utm zone as the projected coordinate system? We inspect pretty small features(10x10ft) and use gps units to mark up these features. I've used both datums and haven't noticed any sort of difference out in the field, so I'm sure it's not a massive difference but I am curious


texts: OpenLayers OpenStreetMap OSM not rendered

texts:
I am trying to render a part of OpenStreetMap offline, previously loaded as map.osm (OSM XML). I'm starting a localhost and loading the xml from disk.

I tried to use similar example (https://openlayers.org/en/latest/examples/vector-osm.html) changing only the source to my file, yet this didn't work out. Data loaded is pretty much similar to data I see in the example.

I tried a bunch of various approaches, and I still can't render the map. I get no errors, and I can't get what I'm missing.

var vectorSource = new VectorSource({
   format: new OSMXML()
 });

 var xml = await fetch('map.osm').then(res => res.text());

 var features = (new OSMXML()).readFeatures(xml);
 // Features are parsed OK
 vectorSource.addFeatures(features);

 var vector = new VectorLayer({
   source: vectorSource,
   // Using styles, I don't post them here (a lot of code)
   style: function (feature) {
     for (var key in styles) {
       var value = feature.get(key);
       if (value !== undefined) {
         for (var regexp in styles[key]) {
           if (new RegExp(regexp).test(value)) {
             return styles[key][regexp];
           }
         }
       }
     }
     return null;
   }
 });
 map = new Map({
   layers: [vector],
   target: document.getElementById('map'),
   view: new View({
     center: [0, 0],
     maxZoom: 0,
     zoom: 0
   })
 });


texts: Which data formats can store Bézier curves?

texts: I've heard that ESRI .gdb can. Are there other databases or file formats that can as well?


texts: Iterate Through All Combinations of Features by Selection in ArcGIS Python Script

texts:
I have question about ArcGIS and python script.
There are following sample code,
So I would like to change Another Field , not OIDFieldName.
Please let me know how to change Another Field.

import arcpy
import itertools

arcpy.env.workspace = r"C:\path\to\workspace"

base_shapefile = "some_shapefile.shp"

oid_field = arcpy.Describe(base_shapefile).OIDFieldName
feature_count = int(arcpy.GetCount_management(base_shapefile).getOutput(0))

for (i, j) in itertools.combinations(range(feature_count), 2):
    arcpy.Select_analysis(base_shapefile , "combination_%i_%i.shp" % (i, j), '"%s" in (%i, %i)' % (oid_field, i, j))

These sample code is combination for 5C2. I want to create all combiantions 5C3 , 5C4 , 5C5.


texts: ET Geotools Split Lines at Vertices Issue

texts: This might seem like a really simple question but for whatever reason I can't figure it out. How do I run the split line tool in ET Geotools on an entire polyline layer? I can do it for a single feature within the layer but I am unsure of how to apply it to the entire layer.


texts:
I am working on creating a project boundary for a specific area. To do this I merged several data sets  including roads,rivers & railroads (polylines) and lakes, civil townships & congressional townships, (polygons converted to polylines) into a single polyline feature class. I then split the polylines at every intersect in the data. I then created a subset of the data using only the nearest features to the project area to create a boundary.

My problem is my data has some issues with overlap. In some areas the civil townships are coincident with bodies of water or roads and the polyline for both types of feature are overlapping (or close to overlapping) in my subset of the data. I would like to clean up the data so there is only 1 line for each section of the boundary and duplicate boundary lines are deleted. after I get a clean polyline feature I was planning and using the feature to polygon tool to create the the final boundary.

Is there a simple method of cleaning up this data? I already have used topology rules to clean up a few thing such as dangles, but I couldn't really find a rule that applied to my situation.

any advice would be appreciated.


texts: PyQGIS standalone script: how to join attributes by location?

texts:
I'm writing a stand-alone script using no GUI.

Earlier I posted a question about having some code issues using JoinAttributes function. However, I found out that the function is not appropriate for my needs. I'm looking for a function that would add attributes of a poligon layer features to point features that lay on it (exactly like the function Join attributes by location in QGIS). I found an earlier answer  suggesting to use modules shapely and fiona. However, I would like to avoid external libraries in this application.

Is there any PyQGIS function that would do the trick in a stand-alone script?


texts: Why is this .prj file causing the shapefile to render in meters instead degrees?

texts:
So I've got a shapefile that has polygons of countries, with coordinates in decimal degrees:

 Country (String) = Algeria
 POLYGON ((6.928073999999953 36.88360652,
           6.949973999999969 36.883306519999991, ...

gdalsrsinfo on the .prj file outputs:

PROJ.4 : '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs '

PROJCS["Mercator",
    GEOGCS["GCS_unnamed ellipse",
        DATUM["unknown",
            SPHEROID["Unknown",6378137,0]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    PARAMETER["latitude_of_origin",0.0]]

To my untrained eye, the GEOGCS[] bit looks like it's defining the coordinate system in degrees, so this should work.  However, if I include the .prj file, Google Maps Engine thinks the coordinates are in meters and renders everything as a tiny blob around (0,0).  If I drop the .prj file, it defaults to degrees and returns a much more sensible rendering.

What's wrong, and how do I fix it?  (Note the the shapefile is from a third party and produced with an unknown toolset.)


texts:
I should move some vertices of a polygon to fix the boundaries of a country, how can I do it while keeping the adjacent polygons still connected? Is there a way to anchor permanently  overlapped vertices of adjacent polygons?

I'm using Qgis


texts: MODIS REST API transform data to GeoTIFF with Python

texts:
I'm trying to transform the response from a request with land surface temperature data from MODIS API to a GeoTIFF:

https://modis.ornl.gov/rst/api/v1/MOD11A2/subset?latitude=0&longitude=0&band=LST_Day_1km&startDate=A2001001&endDate=A2001001&kmAboveBelow=1&kmLeftRight=1

This request returns the following json:

{
  "xllcorner": "-9833349.10",
  "yllcorner": "1575263.27",
  "cellsize": 926.6254330558338,
  "nrows": 3,
  "ncols": 3,
  "band": "LST_Day_1km",
  "units": "Kelvin",
  "scale": "0.02",
  "latitude": 14.182305668542996,
  "longitude": -91.20138157154014,
  "subset": [
    {
      "modis_date": "A2020001",
      "calendar_date": "2020-01-01",
      "band": "LST_Day_1km",
      "tile": "h09v07",
      "proc_date": "2020010222551",
      "data": [15263, 15232, 15219, 15234, 15217, 15219, 15224, 15227, 15299]
    }
  ]
}

I have a function in Python that receives the JSON from before as input and creates a GeoTIFF:

from osgeo import gdal, osr

def createGeoTif(data):
    values = np.array(data['subset'][0]['data'])
    values.shape=(data['nrows'],data['ncols'])
    values = (values * float(data['scale']))-273.15

    drv = gdal.GetDriverByName("GTiff")
    ds = drv.Create('./test.tif',data['ncols'],data['nrows'],1,gdal.GDT_Float32)

    xllcorner = float(data['xllcorner'])
    yllcorner = float(data['yllcorner'])
    x_resolution=y_resolution = data['cellsize']

    ds.SetGeoTransform([xllcorner, x_resolution, 0, yllcorner, 0, y_resolution])

    epsg = 4326

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    dest_wkt = srs.ExportToWkt()

    ds.SetProjection(dest_wkt)

    ds.GetRasterBand(1).WriteArray(values)

The problem is that the GeoTIFF is placed far away where it should be. The problem I think is that the request returns the xllcorner and yllcorner in MODIS sinusoidal projection but I can't find in Google to which EPSG number it belongs so I just need to change this part of the code to fix it:

    epsg = 4326 #Modis sinusoidal code?

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    dest_wkt = srs.ExportToWkt()

    ds.SetProjection(dest_wkt)

From modis api doc:

xllcorner --- MODIS sinusoidal x-coordinate of the lower left corner
of the subset
yllcorner --- MODIS sinusoidal y-coordinate of the lower left corner
of the subset


texts: osm2pgrouting hangs on "classes table created"

texts:
I'm trying to run osm2pgrouting and on Ubuntu 12.04...everything looks like it's working until it gets to "classes table created" when it just stops. I'm watching the process with the system monitor and it looks like it's loaded the whole OSM file into memory(the RAM used by the program is about the same size as the file, ~200MB). But once it gets to that part the CPU usage goes down to nothing and it doesn't do anything. The program doesn't crash or anything, it just seems to go on break.

Here's my command:

osm2pgrouting -file file.osm -dbname osm -user name -conf mapconfig.xml -clean -passwd xxx

Help?


texts:
I have an issue with arcmap. In the layout view and dataview it shows all buildings normal (it's a bit blurry, but it's only for a study not for presentation), but when I export it, all gaps between the buildings get filled. I don't know what the error could be and first thought it would be in some export setting, but it isn't. I made an EMF and inserted the picture in Arcmap and the EMF was filled in too. Then I have been trying to fix or find an answer for 3 hours now, while my Arcmap kept freezing too, but still to no avail. If someone could help me fix this, I would be very grateful :)

Probably there is an easy fix to the problem.

Some more background information:
All files are shapefiles (vector)
When I export just the layer with the building information all works fine.


texts: How to Select By Attributes on only Field with all its chars in Capital letters?

texts:
Need some assistance with this situation: I have a polygon feature class with a string Field. It has mixed values like: "PARIS" - "LISBON" - "london" - "luxemburg" - "Athens".

In this example, I'd only want to SELECT: PARIS and LISBON. How could be the statement? I think it is something with UPPER ("Name_Field") within the Select By Attributes text box... but I don`t know.


texts: Put LeafletJS geoJSON's into a MySQL database through PHP

texts:
I'm doing a project that will let users manipulate a map through LeafletJS, creating circles and polygons on top of the map. I need to gather these inputs for later analysis, using a MySQL database for the storing.

I was thinking in getting the input as geoJSON, then converting it into PHP by separating into several variables to correctly fill the tables in the database.

My question is more about where do I do that conversion... I imagine it is not everything on the client side right? I need to send the geoJSON to some back-end making the JSON-PHP conversion there, right? Or is there any other more direct way to do it?


texts: Using ESRI Dojo version for building web application

texts:
If I am using ESRI's JavaScript API for my application, then I cannot load a different version of Dojo other than the one shipped with ESRI. I have to rely on ESRI for Dojo. The map appears in only 3 pages. Do I have to use ESRI Dojo on the rest of the pages?

Does ESRI dojo load all components or only things specific to ESRI map? If I want to use other components, can I load those from ESRI?


texts: ArcGIS hillshade creates artifact on tile boundary

texts:
I'm working with ArcGIS on a DTM made of 1x1km tiles (cell size 20cm). When I try to create an hillshade model based on this DTM I get visual artifacts on tile junctions, whereas there are no visible drop or gap on the DTM itself.

Edge artifact between tiles

DTM Cross section

I get the same result if I combine first the tiles and create the hillshade model or if I create the hillshade model tile by tile.
I also get the same result if I create the hillshade model with Global Mapper.

Do you have any advice to overcome this issue ?


texts: ArcGIS Python syntax error using nested if else

texts:
I am trying to use the field calculator in ArcMap 10.2 to fill an empty field based on the contents of two exiting fields usint this python script in the Pre-Logic Script Code.

def myfunc(legend,biltype):
    if(legend == '0000 Manmade surface or step'):
        return 'Manmade surface or step'
    elif(legend == '0000 Multiple surface (garden)'):
        return 'Garden'
    elif(legend == '0000 Natural surface'):
        return 'Natural Surface'
    elif(legend == '0000 Path'):
        return 'Path'
    elif(legend == '0000 Railway'):
        return 'Railway'
    elif(legend == '0000 Road'):
        return 'Road'
    elif(legend == '0000 Road traffic calming'):
        return 'Road Traffic Calming'
    elif(legend == '0000 Structure'):
        return 'Structure'
    elif(legend == '0000 Track'):
        return 'Track'
    elif(legend == '0000 Unclassified (or broken)'):
        return 'Unclassified (or broken)'
    elif(legend == '0000 Unknown surface'):
        return 'Unknown Surface'
    elif(legend == '0321 Archway'):
        return 'Archway'
    elif(legend == '0323 Glasshouse'):
        return 'Glasshouse'
    elif(legend == '0377 Boulders'):
        return 'Boulders'
    elif(legend == '0379 Coniferous trees'):
        return 'Coniferous Trees'
    elif(legend == '0380 Coniferous - scattered'):
        return 'Coniferous - scattered'
    elif(legend == '0382 Marsh reeds or saltmarsh'):
        return 'Marsh reeds or saltmarsh'
    elif(legend == '0384 Nonconiferous trees'):
        return 'Nonconiferous trees'
    elif(legend == '0385 Nonconiferous - scattered'):
        return 'Nonconiferous - scattered'
    elif(legend == '0386 Orchard'):
        return 'Orchard'
    elif(legend == '0390 Rough grassland'):
        return 'Rough grassland'
    elif(legend == '0392 Scrub'):
        return 'Scrub'
    elif(legend == '0395 Upper level communication'):
        return 'Upper level communication'
    elif(legend == '0400 Inland water'):
        return 'Inland water'
    elif(legend == '0321 Building'):
        if(biltype.startswith('Commercial, Education'):
            return 'Education'
        elif(biltype.startswith('Commercial, Emergency and Rescue Service'):
            return 'Emergency Services'
        elif(biltype.startswith('Commercial, Retail'):
            return 'Retail'
        elif(biltype.startswith('Residential'):
            return 'Residential'
        elif(biltype.startswith('Commercial, Industrial'):
            return 'Industrial'
        elif(biltype.startswith('Commercial, Hotels'):
            return 'Hotels'
        elif(biltype.startswith('Commercial, Community Services'):
            return 'Community Services'
        elif(biltype.startswith('Commercial, Medical'):
            return 'Medical'
        elif(biltype.startswith('Commercial, Leisure'):
            return 'Leisure'
        elif(biltype.startswith('Commercial, Offices'):
            return 'Offices'
        elif(biltype.startswith('Commercial, Utilities'):
            return 'Utilities'
        else:
            return 'Building - Other'
    else:
        return 'error'

and then setting the field to equal

myfunc( !LEGEND!, !blpu_class!)

Where LEGEND, and blpu_class as the two field with the source data.

When I run the field calculator I get an 000989 Python syntax error:parsing error SyntaxError: invalid syntax (line 51)

This is the start of the nested if, but I can't work out what the problem is.


texts: Feature class to feature class on huge data - iterate loop in ArcGIS pro?

texts:
I have a very large dataset containing a very large area. I have two files; one that is area specific shapefile and a standalone table that I have joined to the shapefile. Earlier, I have run feature class to feature class to make the join permanent but in this instanse, my computer runs out of memory way before I even get close.

What I want to do, is to make ArcGIS Pro run on smaller areas at a time and collect it all in a permanent file. For example; I have county as another feature and I thought I could do something like this:

Run feature class to feature class on the original file join with the standalone table by county 1 and then make the output, then by county 2, then by county 3 and so on. I think I need to make it run over specific selections but I can't seem to make it run automatically over all the counties and then make it in one big file.

Can anyone help me in understanding how to make this work?


texts: GeoServer installation issue on Windows 10

texts: I installed GeoServer 2.13.4 with JRE 1.8.0_241 on Windows 10 64-bit system. The installation was fine however, when I goto http://localhost:8080/geoserver from the browser there is an error 404. I am not sure what the problem is followed the installation steps in documentation for the respective version. How can I resolve the issue?


texts: Extract lat/long from native PostgeSQL point type column

texts:
In my table I have a column location that is a point. I tried this:

SELECT ST_X(location), ST_Y(location) FROM locations;

But I still get the error message:

ERROR: function st_x(point) does not exist


texts: Seeking population of each Indian village and town?

texts:
There are over 600,000 villages in India.
I want to get the population of each of the village according to 2011 census.
I searched the http://censusindia.gov.in website for hours but could not get data for all the villages in a single file (of any format) or a few files.

In the population finder page, I could search for villages one by one, but could not get the data for all the villages in less time.

Can you inform me where to get the data?


texts: Seeking UK Postcode out-code list containing county?

texts:
I’m trying to put together a spreadsheet containing the following fields: postcode out-code, county

Can anybody please point me in the direction of a list (or 2 lists which I could join) containing UK postcode out codes and county?

I need to know for which county a postcode out code is located in.

Preferably open source data.


texts:
I am working on a map to overlay precipitation lines onto Africa, the end goal being to interpret plant growth in different regions. I got a map up, and most of the precipitation lines look fine, but in some cases there are straight lines cutting through portions of the map that probably should not be there.

I am more or less new at this stuff.

Here's the map output:

And here is the code that got me there:

precip <- readOGR(dsn="PrecipSHP", layer = "precip_contours",
stringsAsFactors = FALSE, verbose = FALSE)
precip_f <- fortify(precip)

africa_countries <- readOGR(dsn="Africa", layer= "AfricanCountries",
stringsAsFactors=FALSE, verbose=FALSE)
africa_main <- subset(africa_countries, Land_Type == "Primary land")
country_names <- africa_main@data$COUNTRY
center_pts <- gCentroid(africa_main, byid=TRUE)
labels <- cbind.data.frame(country_names, center_pts@coords[,1],
center_pts@coords[,2])
names(labels) <- c("country","x","y")

africa_f <- fortify(africa_main, region = "COUNTRY")

africa <- get_map(location = c(-19,-36,52,38), zoom = 4, maptype =
"watercolor")
africa <- ggmap(africa)
africa <- africa +
  geom_map(data=africa_f, map=africa_f, aes(x=long, y=lat, map_id=id), color
="white", fill ="orangered4", alpha = .4, size = .2)
africa <- africa + annotate("text", x = labels$x, y = labels$y, label =
labels$country, size = 2) +
  geom_path(data = precip, aes(x = long, y = lat))

This comes out of a larger program, but I tested the code above alone and it worked to produce the map above. You can see that some of the lines, especially in the south, look a bit odd.

Why is that?


texts:
After a LiDAR flight, we can get a file containing the airplane trajectory in SBET format - Smoothed Best Estimate of Trajectory.

Is it possible to open it in a GIS software like QGIS or ArcGIS?

If not, is there a way to convert it into a readable file?


texts: Adding raster layers of different shape using rasterio

texts:
Is there a way to sum two rasters of different shape using rasterio?  To explain it better, the rasters have different bounding boxes, but there is an area that overlaps.

The objective is to get the result of the operation and write it to file, limiting the output raster to the area that overlaps.


texts: Searching for UN regions/sub-regions shapefile

texts:
for a visualisation project I need the boundaries of the UN regions as defined here : http://unstats.un.org/unsd/methods/m49/m49regin.htm

That's 21 polygons for the sub-regions.

I have been looking around on the usual places but couldn't find anything that also includes the oceans. The 21 polygons should cover the whole world so that any point over the world can be associated to one of the regions.

Does someone knows where I could find such file ? Or should I build it... Basically is there an official boundary in the middle of the atlantic ocean that separates the "Western Africa" region form the "South America" region ?

My assumption is that the sub-region contains more than just the countries (for example the "Caribbean" region also contains the ocean between the islands.

To clarify, I know how to identify if a point is in a polygon, that's not the point of this question. I'm looking for a shapefile (or geojson, or anything I can work with) with 21 polygons corresponding to the UN definition and covering the whole world to know if a point located somewhere in the atlantic belongs to west Africa, or to south America, or the Caribbean, etc.


texts: Adding newly created layer to Map in CURRENT project using ArcPy for ArcGIS Pro?

texts:
In this script, I create a feature layer from a file geodatabase feature class and wish to add the layer into a map within the current session in ArcGIS Pro 2.0.  The script crashes on the last line.  How do I fix this error?

aprx = arcpy.mp.ArcGISProject("current")
aprxMap = aprx.listMaps("MainMap")[0]
featureLayer = arcpy.MakeFeatureLayer_management(splitSegmentFC, "SegmentSplit")
aprxMap.addLayer(featureLayer)

The error message is:

Traceback (most recent call last):
File "E:\ProvingGrounds\SCAT\Script\ArcGISpro\ScatTrackProcessing_v10_PRO.py", line 1520, in <module>
aprxMap.addLayer(featureLayer)
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\utils.py", line 191, in fn_
return fn(*args, **kw)
File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 1055, in addLayer
return convertArcObjectToPythonObject(self._arc_object.addLayer(*gp_fixargs((add_layer_or_layerfile, add_position), True)))
ValueError: SegmentSplit
Failed to execute (Script).


texts: Altering sun azimuth on shaded relief LiDAR maps in QGIS?

texts:
I have a LiDAR map and want to change the view so that the ground shadow is manipulated (looking for archaeological anomalies). I think it is the sun angle that I need to alter but I don't know how.

I am running on Windows 7 with QGIS 2.4 downloaded. I am new to QGIS and I am struggling with such task.


texts: Precondition in Modelbuilder - use in next process

texts:
In modelbuilder I am trying to create a tool which will check the coordinate system of the input shapefile. If it is GCS_WGS_1984, then nothing happens. If it is not, then it will reproject it to GCS_WGS_1984.

I know I have to create a script and import the tool to modelbuilder, which I have done. Below is the script I used:

import arcpy, os
from arcpy import env
env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

Footprints = arcpy.GetParameterAsText(0)

desc = arcpy.Describe(Footprints)
type = desc.spatialReference

try:
    if type.name == "GCS_WGS_1984":
        arcpy.SetParameterAsText(1, "True")
        arcpy.SetParameterAsText(2, "False")
    else:
        arcpy.SetParameterAsText(1, "False")
        arcpy.SetParameterAsText(2, "True")

except Exception as e:
AddPrintMessage(e[0], 2)

And my model:

I connected the "NOT" to the "PROJECT" and set it as a precondition, but when I try to run it, I get an error telling me that no input is set for the Project Function. I know that the Footprints are not GCS_WGS_1984. So how do I get the model to run properly?

Is there a problem with the model, script, or both? Thanks!

Ok, so I added in the project tool as suggested, but then I get the following error:

ExecuteError: ERROR 000622: Failed to execute (Project). Parameters are not valid.
ERROR 000628: Cannot set input into parameter out_coor_system.

My script now looks like

import arcpy, os
from arcpy import env
env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

Footprints = arcpy.GetParameterAsText(0)
NewProject = arcpy.GetParameterAsText(2)
OutCoorSys = "GCS_WGS_1984"

desc = arcpy.Describe(Footprints)
type = desc.spatialReference

if type.name == "GCS_WGS_1984":
    arcpy.SetParameterAsText(1, "True")
else:
    arcpy.Project_management(Footprints, NewProject, OutCoorSys)


texts: FeatureRequest different from different Geoserver instances

texts:
I have two instances of geoserver on different systems, with the same Data underlying. One is version 2.4.4, the other one 2.6.2.

Both behave different when I make a Feature Request:

the first one is giving me Bounding Box Data like this:

<gml:boundedBy><gml:Envelope srsDimension="2" srsName="urn:ogc:def:crs:EPSG::31466">
<gml:lowerCorner>5770863.855 363273.177000001</gml:lowerCorner>
<gml:upperCorner>5770864.5512810005 363274.2701699999</gml:upperCorner>
</gml:Envelope></gml:boundedBy>

The second one provides location Data in a different way:

<gml:Polygon ...>
  .....
  <gml:posList>
        5770863.855  xxx yyy 363273.177000001 5770864.5512810005 zzz qqq 63274.2701699999 rrr sss
  </gml:posList>
 ....
</gml:Polygon>

I want to have it on both instances like in the first one, with the bounded Box.


texts: Can you use GetLegendGraphic for a ArcGIS Server 9.3 WMS?

texts:
The following request should return a legend graphic for a particular layer:

http://myserver/arcgis/services/default/MapServer/WMSServer?Service=WMS&Request=GetLegendGraphic&Version=1.1.1&Layer=3&FORMAT=image/png

However instead I get the following error:

<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://schemas.opengis.net/wms/1.1.1/exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1">
  <ServiceException code="InvalidFormat">
Can't parse XML request.
  </ServiceException>
</ServiceExceptionReport>

The GetCapabilities request for the service for the layer contains the following:

<Style>
<Name>default</Name>
<Title>Layer Name</Title>
<LegendURL width="100" height="62">
<Format>image/png</Format>
<OnlineResource xlink:href="http://server/arcgisoutput/wms/default1.png" xlink:type="simple"/>
</LegendURL>
</Style>

The OnlineResource href is a valid PNG file that is accessible externally.

It seems others have had issues, and recommend using external capabilities files, but when the legend already has a valid image file is this necessary? Has anyone successfully configured a WMS to support the GetLegendGraphic request?


texts: Analyzing hot spots using point attributes?

texts:
I'm interested in looking at where hot and cold spots of certain social indicators (such as poverty, housing characteristics, urban sprawl, etc...) are in coastal communities along the East Coast. I have point and polygon features containing the same attribute information. My point data are centroids of the polygons.
My first thought was to use the Hot Spot Analysis tool in arcgis but, because I'm looking at coastal communities, my polygons don't have many neighbors (some don't have any) so I don't think my results would be very accurate. So I thought there must be a hot spot analysis that  I can do using my centroid points that have the same attribute information but the only analyses I can find have to do with spatial clustering of points- not of their attributes. Does anyone know of an analysis I can do to solve my problem?


texts: How can I add data from points to a raster file?

texts:
I have a raster layer and a shape layer consisting of points. I want to add the data from my points to the raster layer as I need this format.

Is there a way in QGIS (preferably), Python or R to add data to my raster file?
EDIT: Here is the data in my shapefile


texts: Interpolating heights in DEM using ArcGIS API's in C#

texts:
I am new to ArcGIS and not familiar with it. I want to interpolate heights in a DEM file using ArcGIS? Is it possible?
Actually, I want to write a C# code that will read the height from DEM file for given coordinates, and if the height is missing then I will perform interpolation and return the height.
I want to use ArcGIS interpolation algorithms instead of writing my own.


texts: Why when I draw autocomplete polygons using a frame`s border, it doesn`t work inside the frame?

texts:
I try to create an Auto-complete polygon within a bigger frame (from the same layer), using frame's border. The thing is that when I try to create it outside the frame, the polygon appear but when I'm trying to do inside, it doesn't work. I tried many times and every time, outside it worked and outside,no.

Why is this situation ? What can I do in order to succeed making autocomplete polygons inside the frame?


texts:
This python code works in ArcGIS Pro but Not in ArcMap . I read some issues that mention a get around by making a feature layer. I tried that and I then get an error saying polytest doesn't exist. maybe since I am writing to a temp file?

    import arcpy, os, shutil, sys

# Path to store geoprocessing files
user_profile_path = os.environ['USERPROFILE']
gdb = r"{}/AppData/Local/Temp/PolyChecker/PolyCheck.gdb".format(user_profile_path)

# GDB already exists, delete it
if arcpy.Exists(gdb):
    shutil.rmtree(r"{}/AppData/Local/Temp/PolyChecker".format(user_profile_path))

# Create fresh GDB path
os.makedirs(os.path.dirname(gdb))

# Create GDB
arcpy.CreateFileGDB_management(os.path.dirname(gdb), os.path.basename(gdb))

# Script arguments
MapUnitPolys = arcpy.GetParameterAsText(0)
if MapUnitPolys == '#' or not MapUnitPolys:
    MapUnitPolys = "MapUnitPolys"  # provide a default value if unspecified

MapUnitPolys_CopyFeatures = arcpy.GetParameterAsText(1)

# Set Geoprocessing environments
MapUnitPolys = MapUnitPolys

# Validate that all Polygons have a map unit

invalid_polygon_found = False

with arcpy.da.SearchCursor(MapUnitPolys, ['SHAPE@', 'MapUnit', 'OBJECTID']) as cursor:
    for row in cursor:
        #arcpy.AddMessage(str(row[1]))
        # Does this Polygon have a map unit
        if row[1] == "" or row[1] == "<Null>" or row[1] is None or row[1] is 0:
            invalid_polygon_found = True
            arcpy.AddMessage('Polygon OBJECT ID:{} is missing map unit... exiting.'.format(row[2]))

# Invalid polygons were found, terminate
if (invalid_polygon_found):
    sys.exit(1)

Polygon_Neighbors = "{}/polytest".format(gdb)

PolygonNeighbor_TableSelect = "{}/PolygonNeighbor_TableSelect".format(gdb)

inFeatures_lyr = "{}/inFeatures_1yr".format(gdb)

# Process: Polygon Neighbors
arcpy.PolygonNeighbors_analysis(MapUnitPolys, Polygon_Neighbors, "OBJECTID;MapUnit", "NO_AREA_OVERLAP", "BOTH_SIDES",
                                "", "METERS", "SQUARE_METERS")

# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management(Polygon_Neighbors, "NEW_SELECTION", "src_MapUnit = nbr_MapUnit")

# Process: Table Select
arcpy.TableSelect_analysis(Polygon_Neighbors, PolygonNeighbor_TableSelect, "src_MapUnit = nbr_MapUnit")

arcpy.GetCount_management(PolygonNeighbor_TableSelect)

arcpy.AddMessage(arcpy.GetMessages())

if int(arcpy.GetCount_management(PolygonNeighbor_TableSelect)[0]) > 0:
    arcpy.MakeFeatureLayer_management(MapUnitPolys, inFeatures_lyr)
else:
    print ("done")

# Process: Add Join
arcpy.AddJoin_management(inFeatures_lyr, "OBJECTID", PolygonNeighbor_TableSelect, "src_OBJECTID", "KEEP_COMMON")

# Process: Copy Features
arcpy.CopyFeatures_management(inFeatures_lyr, MapUnitPolys_CopyFeatures, "", "0", "0", "0")

# Process: Remove Join
arcpy.RemoveJoin_management(inFeatures_lyr, "")

# Execute Delete
arcpy.Delete_management(PolygonNeighbor_TableSelect)
arcpy.Delete_management(Polygon_Neighbors)

arcpy.AddMessage('All done! Check Polygons')

error1

    Traceback (most recent call last):
  File "C:\Downloads\GeMS-ArcPro-Toolbox-master\GeMS-ArcPro-Toolbox-master\Scripts\PolyChecker2.py", line 56, in <module>
    arcpy.SelectLayerByAttribute_management(Polygon_Neighbors, "NEW_SELECTION", "src_MapUnit = nbr_MapUnit")
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\management.py", line 7726, in SelectLayerByAttribute
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
The value cannot be a table
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (SelectLayerByAttribute).

Failed to execute (PolyChecker).
Failed at Mon Aug 12 13:28:48 2019 (Elapsed Time: 3.90 seconds)


texts:
With the following code I am trying to produce a NDVI graph with the days of the year on the x axis (based on this example, 2nd graph).

var modisNDVI = ee.ImageCollection ('MODIS/MCD43A4_006_NDVI');
    var modiscollection = ee.ImageCollection (modisNDVI.filterDate('2007-05-01','2011-09-30'));
    var clip=modiscollection.mean().clip(region);
    var chart = ui.Chart.image.doySeriesByYear(
        clip, 'NDVI', region, ee.Reducer.mean(), 500);

    print(chart);

    Map.addLayer (clip, {min:0.0,max:1,palette:['FFFFFF','CC9966','CC9900','996600','33CC00','009900','006600','000000']},'NDVI');

From this, I get the following error and as a new user, I can't understand it.

P.S. For the region I just draw a polygon

Error :

Collection.first: Error in map(ID=0):
    Date: Parameter 'value' is required.


texts: Limiting number of decimal digits in ArcGIS to have exported file in CAD with same format?

texts: I am supposed to have some numeric values in ARCGIS 10.2 with 3 decimals Digits and then convert them to Microstation or AutoCAD. I did this Setting using field calculator and 'round' function in arcgis but when I convert it to CAD again it would come with more than 10 decimals which is not desired.


texts: I want to make simple regular point grid for a polygon but of known distance in meters between points, in QGIS using regular points research tool. How do I input my own distance? The spacing value '0.0001' is not clear what metric that is


texts:
QGIS 2.18.5

The radius and positioning of my viewsheds seem to be doing something peculiar. I have set my CRS to British National Grid for both the terrain raster (downloaded from USGS) and my observation points. On the fly projection is unchecked. From the image below you see I have 2 observation points, one by the coast to the south, another to the north further inland.

My settings in the viewshed analysis plugin:

The results show both viewsheds out of place. The plugin thinks the observation points are some distance to the south of where they actually are. Not only that, but why is the radius of the calculated area egg shaped? Both of these questions I thought related to a CRS issue. I tried saving the layers to World Geodetic System 1984 and WGS 84. Both attempts led to a python error (the same CRS applied to each layer at the time).
Looking at the uppermost viewshed, why isn't the entire valley shaded in? If I was standing at that point, I'm pretty sure I could see the ground around me...the plugin seems to be underestimating the results, or perhaps I'm using the plugin incorrectly.


texts:
I have a search cursor that using a SQL query to match addresses and then returns the associated roll number with that address. However, some addresses do not match the roll number and then I believe the search cursor stops. Is there a bit of code I can throw in that will continue with the next address even if a match isn't found? I think an 'if' statement would work well, but I don't know what to write for it!

I'm having trouble generating this.
My code is iterating through foldernames to get the address and querying a parcel layer.

I can post the code later, if you need an example.

import arcpy
import os
import re
import sys
import traceback
import collections

basedir = r"C:\Test"
fc = r"G:\ParcelsFixNames.shp"
field1 = "ADD1"
field2 = "STREETNAME"

try:
#find all foldernames in directory
    for fn in os.listdir(basedir):
        string = str(fn)
#split name at hyphen
        if re.findall('-', fn):
#right side of hyphen is street address, leftside is street name
            lhs, rhs =  string.split('-')
            street = str.lstrip(lhs)
            street = str.replace(street,'.','')
            street = str.rstrip(street)
            street = str.upper(street)
            add = str.lstrip(rhs)
            expression3 = '"{field1}" = \'{add}\' AND "{field2}" LIKE \'%{street}%\''.format(field1=field1, add=add, field2=field2, street=street)
            rows = arcpy.SearchCursor(fc, fields="STREETNAME; ARN; ADD1", where_clause=expression3)
            for row in rows:
                 Roll = row.getValue("ARN")
                 newname = str(Roll)
                 newpath = os.path.join(basedir, newname)
                 #if the name exists already add _1 to end
                 if os.path.exists(newpath):
                    os.rename(os.path.join(basedir, fn),
                        os.path.join(basedir, newname + "_" + "1"))
                 else:
                    os.rename(os.path.join(basedir, fn),
                        os.path.join(basedir, newname))
        else:
            print add + " " + street + " " + "Already exists...passing"
            continue
except:
    print ("Error occurred")
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n    " + \
        str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"
    msgs = "GP ERRORS:\n" + arcpy.GetMessages(2 )+ "\n"
    print (pymsg)
    print (msgs)


texts: How do I edit line and fill setting when using Add Feature?

texts:
Newbie here. I see that the creation of new shape layers must be possible, but I cannot work out how to do so without QGIS only using the default ugly green fill polygon setting. Is there a way to create non-ugly polygons in QGIS please, ideally with no fill so I can actually see what is underneath?

I have searched but all the answers appear to be in Python code, rather than "you do it like this" answers.


texts:
I can't show these icons from FontAwesome:

var map = new ol.Map({
        layers: [
          new ol.layer.Tile({
                source: new ol.source.OSM()
            })
        ]
        , target: 'map'
        , controls: ol.control.defaults({
            attributionOptions: {
                collapsible: true
            }
        }), //Vista inicial del mapa
        view: new ol.View({
            center: ol.proj.fromLonLat([-5.9334, 43.5848])
            , zoom: 15
        })
    });
    var vectorLayer = new ol.layer.Vector({
        source: new ol.source.Vector({
            features: new ol.format.GeoJSON().readFeatures(geojsonObject, {
                featureProjection: 'EPSG:3857'
            })
        })
        , style: new ol.style.Style({
            /*image: new ol.style.Circle( /** @type {olx.style.IconOptions} ({
                radius: 20
                , fill: new ol.style.Fill({
                    color: '#000000'
                })
            }))*/
            text: new ol.style.Text({
                text: '\uf041'
                , font: 'normal 38px FontAwesome'
                , fill: new ol.style.Fill({
                    color: '#ff0000'
                 })
            })
        })
    });
    map.addLayer(vectorLayer);

Does anybody use typographic icons as Openlayers icons?
I need help...


texts: Determining projection of Albanian military map?

texts:
I've received some very useful maps of Albania for my work (scale 1:25000) and I would like to georeference them in order to put them together to one Map for my destined regions.

Unfortunately I cannot find the projection of the maps, they don't seem to be written on the map (it is unfortunately in Albanian language).

What projection might it be?

I'm using QGIS.

full image


texts: I have .dwg files which are not georeferenced. Is there any way to add those to Cesium Globe and extrude?


texts: Testing species distribution model residuals for spatial autocorrelation?

texts:
I build a species distribution model with data from location A and tested it on location B. The picture shows a rastermap of location B and what I guess should be the model residuals.

Green - are the areas where the model was right (value 0)
Blue - are the areas where the model predicted species occurrence but there was no observation (value -1)
Black -  are the areas where the species was observed but the model predicted no occurrence (value 1)

I think the clustering is obvious but I also want to measure the spatial auto correlation and give it a value. In Arcgis and R I found the morans I functions. However both need  point data as input.

Should I transform the Raster to points (one point each cel) and use this as basis for the calculation?

What should I do with the areas where the model was right (value 0)?

should I delete them before I perform the test?

This question is related to Calculating residuals of species distribution model and test for spatial auto correlation at Cross Validated.

Edit:

I asked the question because I always read that you have to take care of spatial autocorrelation in the analysis of species distributional data when building a SDM:

"Species distributional or trait data [...] often display spatial autocorrelation, i.e. locations close to each other exhibit more similar values than those further apart. If this pattern remains present in the residuals of a statistical model based on such data, one of the key assumptions of standard statistical analyses, that residuals are independent and identically distributed (i.i.d), is violated. "(Dormann, 2007).

However I cannot figure out whether it is an issue in my case. I did an systematic, full sampling of species occurrence/absence and environmental variables in Area A. This data (not georeferenced) I used to build a deterministic model to predict species distribution based on habitat preference. I just built habitat preference curves for each predictor based on species occurrence and habitat availability at each sample point and combined the habitat preferences for each predictor to a combined preference by multiplication. @Jeffrey Evans: this could also be interpreted as occurrence probability.

Then I applied the model to the independent area B, where I have a grid/raster of the environmental variables. I also mapped the species occurrence in B to test the model against independent data. I tested the occurrence probability against the observed occurrences. My AUC values are good (around 0.8).

Now I wanted to see whether there is some kind of spatial pattern in the residuals of the model. Iam not even sure whether I should call it residuals because i didn't build a statistical model and do the test against independent data. Maybe better call it error, the error between observed and predicted occurrence. My understanding is that when there is an error which shows spatial autocorrelation then my model omits an important predictor or process that explains the spatial distribution

My main worry is whether I can I trust the AUC value and assess the predictive performance as good?

Because I could not figure out on which data I should do the morans I test I askesd the initial question.


texts: Line between outside and inside polygon - geotools


texts: Two boundary sets intersect - some are genuine diffferences, others slivers resulting from poor digitising. How to differentiate?

texts:
Problem:

I have two boundary sets. Intersecting reveals that there are genuine differences between them - as expected. However, there are also a host of slivers due to their being digitised by different agencies. These mismatches may have statistical implications.

I have tried several approaches to clearly differentiate one mismatch type from the other.

Current thinking:

I get the MBR of each polygon in boundary set 1 and get the ratio of the area of the MBR to the area of the actual polygon. I do the same for boundary set 2.
I get the MBR of each intersection (these have been disaggregated to singlepart) and get the ratio of the area of each intersection's MBR to the area of the actual intersection.
The distribution of the ratio for slivers is likely to be materially different to the distribution of the ratio for genuine differences.

QUESTIONS

Do you agree that this seems like a sensible approach?

How easy is it to obtain the MBR (minXmaxXminYmaxY) for a polygon in ArcGIS 10.1 Advanced? Either point-and-click or Python suggestions welcome.

If my thinking is flawed, please point out how, and suggest an alternative.


texts: Proportional Symbols - Custom legend breaks

texts:
I mapping the number of Boko Haram victims in each state of Nigeria for several years. I chose to make a proportional symbol maps using the centroid of each state. For a given year, the area of the circular symbol is directly related to the number of deaths.

Using ArcMap's proportional symbols symbology does the job, but I have a problem with the legend:

The default legend offers me 3 categories which are not very convenient (1, 100, 10.000). I would like to be able to change the breaks in the legend to, say, 100, 1.000, 5.000, 10.000.

Is there a way to do this?


texts: QGIS - is it possible to convert accented characters to regular ones in attribute table?

texts:
I have OSM dataset containing names with national characters (this is CSV export, but I have it in SHP):

osm_id,code,fclass,population,name,longitude,latitude
276645544,1003,village,221,Krč,14.248410,49.194988
444571400,1003,village,88,Radonice,14.548085,49.143251
488567474,1003,village,115,Doubravka,14.363835,49.282158
488567478,1003,village,114,Koloměřice,14.399437,49.280328
488567484,1003,village,111,Údraž,14.309249,49.281551
491452727,1004,hamlet,33,Poněšice,14.481200,49.105121
601855465,1003,village,147,Hvožďany,14.433230,49.286094
631878480,1003,village,211,Myšenec,14.208999,49.223213

and would like to converte them to this form:

osm_id,code,fclass,population,name,longitude,latitude
276645544,1003,village,221,Krc,14.248410,49.194988
444571400,1003,village,88,Radonice,14.548085,49.143251
488567474,1003,village,115,Doubravka,14.363835,49.282158
488567478,1003,village,114,Kolomerice,14.399437,49.280328
488567484,1003,village,111,Udraz,14.309249,49.281551
491452727,1004,hamlet,33,Ponesice,14.481200,49.105121
601855465,1003,village,147,Hvozdany,14.433230,49.286094
631878480,1003,village,211,Mysenec,14.208999,49.223213

where accents are removed so instead of "č,š,ř,á" etc. you have "c,s,r,a". I does not have to be editing of current column - could be creating new column with converted names.

I am currently going to use some online tool like this one for the CSV file but it would be great if I could solve this in QGIS directly.


texts: How to Calculate Utilization Distribution Overlap Index (UDOI) in ArcGIS Desktop?

texts:
I was wondering if any of you happen to have any information on how to calculate the Utilization Distribution Overlap Index (UDOI) in ArcGIS 10? If I understand it correctly the UDOI takes two kernel density maps and then calculates the percent overlap between the two using a range between 0 and 1. 1 indicating 100% overlap while 0 indicates no overlap. For those that have access here is the journal article outlining the concepts behind UDOI.

http://onlinelibrary.wiley.com/doi/10.2193/0022-541X(2005)69[1346:QHOTIO]2.0.CO;2/abstract

In the literature UDOI is always calculated using the adehabitat package in R. However, I have zero R knowledge and the code is a bit daunting to say the least. It seems to me that this could be done using tools that are already available in Arc.


texts:
When executing: suite-sdk deploy procedure I keep getting the following errors:

Checking provided application path.
Building application resources for 'C:\xxx\xxx\Sites\mySite'.
Aug 04, 2015 1:08:45 PM com.google.javascript.jscomp.LoggerErrorManager println
WARNING: input.js:79235: WARNING - Parse error. Non-JSDoc comment has annotations. Did you mean to start it with '/**'?
    /* @preserve

Aug 04, 2015 1:08:45 PM com.google.javascript.jscomp.LoggerErrorManager println
SEVERE: input.js:72815: ERROR - Parse error. Internet Explorer has a non-standard intepretation of trailing commas. Arrays will have the wrong length and objects will not parse at all.
        }],
         ^
...few more of the same errors...

Aug 04, 2015 1:08:45 PM com.google.javascript.jscomp.LoggerErrorManager printSummary
   WARNING: 6 error(s), 1 warning(s)
    Error: Unable to compile source. (C:\Program Files (x86)\OpenGeo\sdk\lib\buildkit\lib\buildkit\closure.js#110)
at C:\Program Files (x86)\OpenGeo\sdk\lib\buildkit\lib\buildkit\closure.js:110 (anonymous)
at C:\Program Files (x86)\OpenGeo\sdk\lib\buildkit\lib\buildkit\closure.js:137 (anonymous)
at C:\Program Files (x86)\OpenGeo\sdk\lib\buildkit\lib\buildkit\build.js:56 (main)
at C:\Program Files (x86)\OpenGeo\sdk\lib\buildkit\lib\buildkit\build.js:73

BUILD FAILED
C:\Program Files (x86)\OpenGeo\sdk\build.xml:209: Java returned: -1

Total time: 8 seconds

Does anyone have an idea what is causing this? I  cannot find any decent answers on Google regarding this!


texts: Getting from wkt to shapefile in Python?

texts:
I am querying tables from a SQL-Server database and would like to convert the table into a shapefile.  It has wkt geometries and wkb as columns in the table, along with other attributes that need to be fields in the dbf.  I'm wondering what the most efficient way to convert these returned results into shapefiles is.  I'm not sure ogr2ogr is the best method because I already store the attributes in a list of dictionaries in my script (for sorting and filtering purposes).  I've been looking at pysal and fiona, but can't find any examples showing exactly what I'm trying to do, and I'm rather new to Python.  Also, I do have access to arcpy.

Update after @fenris comments:

I've tried this method but am still getting errors.

Here's what I'm doing:

from shapely.geometry import mapping
from shapely.wkb import loads
from fiona import collection
outDir = r'C:/Users/x/Documents/'
outFile = 'Output_CIM'
outCellShp = outDir + outFile + '.shp'

schema = {'geometry': 'Point', 'properties': {'veh' : 'str',
               'id' : 'str',
               'imStart' : 'datetime.datetime',
               'imEnd' : 'datetime.datetime',
               'area' : 'Decimal',
               'c' : 'Decimal'} }

with collection(outCellShp, "w", "ESRI ShapeFile", schema) as output:
    print('here')
    for row in images:
        geometry = loads(row['wkb'])
        output.write({'properties': {'veh' : row['veh'],
               'catid' : row['id'],
               'imStart' : row['imStart'],
               'imEnd' : row['imEnd'],
               'area' : row['area'],
               'cc' : row['c']},
               'geometry' : mapping(geometry)})

And am getting the following error:
Error: class 'fiona.errors.DriverError'

It doesn't appear that it's making it passed the collection() call.  Am I allowed to use 'datetime.datetime' as a properties type? I'd actually prefer to have that as a string because Arc cuts it off at day anyway.  Right now I'm using wkb as a bytearray, but also have AIS_Geometry and AIS_Geography in the database, if either of these would yield better results.

I was able to get the above working with the following:

schema = {'geometry': 'Polygon', 'properties': {'veh' : 'str',
               'id' : 'str',
               'imStart' : 'str',
               'imEnd' : 'str',
               'area' : 'float',
               'c' : 'float'} }

with collection("outCellShp.shp", "w", "ESRI Shapefile", schema) as output:
    for row in images:
        geometry = loads(bytes(row['wkb']))
        output.write({'properties': {'veh' : row['veh'],
               'id' : row['id'],
               'imStart' : str(row['imStart']),
               'imEnd' : str(row['imEnd']),
               'area' : float(row['area']),
               'c' : float(row['c']}),
               'geometry' : mapping(geometry)})

I had to make sure there was a polygon in the schema as opposed to a point.  Additionally, I had to convert the bytearray into a bytestring, the decimals into floats, and the datetime.datetime into str.


texts: Importing ESRI shp / gdb files with M value into Postgis

texts:
Has anyone done this successfully? I have been using DB Manager in QGIS 2.12 and my MutilinestringM values are converted to Multilinestring in PostGIS. When I use the PostGIS shapefile importer, this does not happen and the lines are converted successfully.

I thought based on reading this that perhaps better support for M values was coming in GDAL 2.0. However, I have tested with a nightly release of QGIS using this version and that doesn't seem to be the case.

There are a number of different components to this so it's difficult to know which particular element is at fault here, the plugin, QGIS in general, GDAL or even PostGIS.


texts:
I have a map with two main components:

a polygon containing a number of shapefiles; and
a colored basemap (from ArcGIS) that covers the entirety of the polygon and beyond.

Is there a way to "gray out" the portions of the basemap that are outside the boundary of the polygon and retain the colors of the basemap within the polygon area?

My end goal is to highlight and draw the viewer's attention to the polygon.


texts: How to automate flooded area calculation and make report using PyQgis?

texts:
I have 3 shapefiles with me. They are a flooded layer, District layer and Village layer. To estimate flooding village wise, I have to do intersecting the features and get the layer. Then I calculate the flooded layer.

Raw village shapefile contains the attributes 'Village_Name' and 'Area'.

District layer contains 'District_name' and 'Revenue_Circle'

After intersecting the area for the flooded layer is calculated as 'flooded_area'. So the intersected layer will have 'District_nam','Revenue_Circle', 'Village_area' and 'flooded_Area'.

After doing all these I have to convert all the attributes to excel file and using Pivot table all the data need to be sorted district wise. In the excel I have to remove underscore, change the field name, adding a percentage of the village area got flooded and display summation at the last and all these are done in separate sheets for each district.

Is there any automatic workflow for this if I have to do this every day?


texts:
I am trying to follow the example to create a new REST service in Geoserver:

http://docs.geoserver.org/stable/en/developer/programming-guide/rest-services/implementing.html

After building, accessing the service at http://localhost:8080/geoserver/rest/hello.txt doesn't do anything.

(REST is working since http://localhost:8080/geoserver/rest  will show me the Geoserver Configuration API)

I have also tried within eclipse.  To me it seems that I need to add some configuration somewhere to get the rest services to know about my new "hello world" rest service and add it to the geoserver rest services?

Can someone explain how to get this example to work?

Thanks


texts: Mapproxy skip generation of blank images

texts: How to configure mapproxy-seed so that it would not create bank images for areas where there are no image data?


texts: Getting points along line at defined distance using PostGIS?

texts:
I was reading a lot of similar questions here but no one ask about PostGIS (or I didn't find it). I have a line in PostGIS and I want to get the point along the line at defined distance.

Is it possible?

I found st_locate_point but it doesn't work for I want to do.


texts: Join Attribute by location (Summary) not capturing all inputs

texts:
I created a 0.1 acre hexagonal grid and am trying to produce a heat map based on the respective polygon area that are found within a hexagon. I am using the Join Attribute by location (Summary) tool, but am having some problems with the output. Image 1 shows my hexagonal grid. Image 2 shows the inputs for the Join Attribute by location (Summary) tool. Image 3 shows the outputs. I am a little confused here because it is saying that it has skipped a couple features, but I made sure that there were no polygons outside on my hexagonal grid (I performed a clip earlier on). Image 4 shows the output. Image 5 shows the graduated symbol style that I picked and image 6 shows that output. You can see some of the hexagons are disappearing. I believe this is because that hexagon has no value for the summed area of the polygons (NULL). I guess I could manually go in and change all NULLS to zero?? Finally, image 7 and image 8 show a hexagon that does contain a polygon, but is not captured when running the Join Attribute by location (Summary) tool. Is there something that I am missing? Could this be one of the "missing features". What suggestions are there to fix that?
image 1

image 2

image 3

image 4

image 5

image 6

image 7

image 8


texts: Calculating flowaccumulation in specific areas

texts:
I there a way of calculating only the flowaccumulation which are generated in the yellow area and assigning this calculated value to the related end of the stream which is entering the river?

red cells represent my linear stream network
blue cells represent the river
black arrows mark the entering points

I´m using ArcGIS 10.2.


texts:
is there any way to merge/union curved polygons in Postgis?

"SELECT ST_UNION(the_geom) from public.polygons"  can be executed, if the type of the_geom is Polygon, but doesn't work with the type CurvePolygon ("ERROR: Unknown geometry type: CurvePolygon").


texts: Modelbuilder - Join - Keep matching records only - export to feature class

texts: I need to join a table to a feature class in a model.  Have it keep only the matching records and then export that joined feature class into a new feature class.


texts: Improve Catchment polygon processing result accuracy?

texts:
I use the ArcHydro Toolbox and ArcGIS 10.4

I ran through the regular steps of Filling Sinks of a DEM, Flow Direction, Flow Accumulation, Stream Definition, Stream segmentation, catchment grid delineation and finally catchment polygon processing.
My result looks like the following picture.
As you can see the catchments are always a little off. From my understanding they should always cut the segments of a stream nicely. Any idea what might cause this error and what to do about it?


texts:
I'm trying to configure Positions Desktop AddIn (v10.5.0.2) with Esri ArcMap(10.4.1) for the first time yet not having any luck. The Positions single use license is activated and the software is showing as installed on my computer.

In ArcMap I navigated to Customize-> Add-In Manager->Customize->Add from File... then select the Positions AddIn from my hard drive folder location.

Once doing so I get the error message prompted from ArcMap(See attached Picture). Does anyone know if I'm doing some wrong with the configuration for this AddIn? How do I resolve this issue? I couldn't find any literature from Trimble stating compatibility issues with (v10.5.0.2)


texts: What is the equivalent of city, county and state in the Philippines?

texts:
I am developing a map application using HERE Maps and wondering if what is the equivalent of the fields 'city', 'county' and 'state' for the addresses in the Philippines?

Currently, I have mapped these fields as follows:

city   == municipality
county == province
state  == region

(I am not so sure if this question is valid here. :) )


texts: Using addbasemap command in ArcObjects.NET?

texts:
Here I executed the save edit command programmatically

      pUid.Value = "{59D2AFD2-9EA2-11D1-9165-0080C718DF97}"   'save edits
      pCmdItem = pApp.Document.CommandBars.Find(pUid)
      pCmdItem.Execute()

How do I execute the addbasemap command programmatically in ArcObjects.NET?


texts:
Let me explain what I mean by "Deep Merge" and what I am trying to do.

Suppose that I have 3 layers (Layer1, Layer2 and Layer3) with the same field names (field_a, field_b and field_c). For any record, each layer has a value only in one of the fields. For example, one row in Layer1 has a non-empty value only in field_a. One row in Layer2 only in field_b. One row in Layer3 only in field_c. The tricky part is that, different rows in the same layer may not have a non-empty value in the same field. My goal is to merge all 3 layers into one target layer, which contains all non-empty values from all 3 layers.

In my mind, it involves some serious loops or joins. So I am wondering if anybody knows any existing tool (open source or commercial) to do that.


texts:
I have created a multipoint shapefile file and added lat/long data for a few hundred points (by copying and pasting columns of data from excel). The new shapefile is the same coordinate system as my data frame and all the other shapefiles.

How do I make these points display? My attribute table has everything it needs, why is nothing happening?

I tried importing my data from an Excel table but it kept freezing my computer so I'm trying a different method. I'm using ArcGIS 10.2.


texts: I have already loaded osm data into postgresql earlier. Now I want to update the database. I have downloaded the pbf file from geofabrik. How I can update the database without deleting the existing Tables?


texts: Change default data source encoding in QGIS to utf8

texts: Working in languages with special characters, I always struggle with remembering to change the data source encoding from system to utf8 when I import files. Is there a way to have it automatically use utf8 instead of system, or maybe changing the data source encoding of system to utf8?


texts: Proximity query in Hive using ESRI geospatial library

texts: I have one table of lines and one of points and need to find all pairs of line + point with a proximity < 10 meters.  The documentation for the ESRI hive library completely elides any discussion of units for distance, but I'm assuming it to be in fractional degrees in a WGS84 CRS.  Is there any way to get to linear distance using the Hive ST_* functions in an SQL query?


texts:
I have problems to use the function v.what.vect and I'm not sure what is the issue. I have a simple forest polygon and try to determine if a point is in the forest or not.
# check the output of the feature nature
v.db.select nature col=landuse where="cat<5"
# function returns
landuse
forest
forest
forest
forest
# add column to points
v.db.addcolumn map=points col="landuse varchar(50)"
# use function v.what.vect
v.what.vect map=points column=landuse query_map=nature query_column=landuse
# function returns
Finding nearest features ...
 100%
 100%
Update vector attribute...
 100%
2696 categories - no nearest feature found
v.distance complete. 2696 records updated.

For me its looks like everything worked fine, just the that the result 0 is wired. Any idea what went wrong or how to check it?

EDIT
I checked again the notes of the grass manual:

In case of a multipoint input map map, with several points having the same category number, it can happen, that the query result is NULL, if the same category number falls into different query_map polygons.

But that is also given:
v.db.select points col=cat where="cat<5"
# function output:
cat
1
2
3
4
5


texts:
I have imported a 508GB DB with unknown SRID with ogr2ogr, hence my geometry SRID is 900914.

Now I'd like to run

SELECT UpdateGeometrySRID('contours','wkb_geometry',3857);

However, I'm running low on disk space (~300GB left) and it seems that UpdateGeometrySRID consume everything before failing.

Is there a way to run things more progressively ?

Edit my own question with best (easiest) answer:

Re-run ogr2ogr with -a_srs epsg:3857 option and remember it for the next time.


texts: ArcGIS desktop 10.3: Is it possible to have all popups showing at once when exporting/printing?

texts: I'm very new to ArcGIS software, and I need a map to look like this map I drew on in MS Paint, with multiple popups showing.  Is there any easy, user-friendly way of going about doing that with ArcGIS 10.3?  I need to able to print the map with all the popups displayed.


texts:
In which units do I have to calculate Slope to use as an input for the LS Factor tool in SAGA (degrees, %, radians)?

I couldn't find anything in the documentation.


texts: Offsetting 100 ft from line feature while editing in ArcMap?

texts: I have the ArcGIS 10.6 advanced lisense. I am trying to offset a line feature using the editor tool. I know that I can do it in COGO, but how to set the unit as ft. In COGO, the offset unit option is not there.


texts: Add MapBox sourced hillshade to MapBox Studio project

texts:
I am working on a MapBox studio project that was originally the pre-styled 'OSM Bright 2'.

I have since made lots of edits and would now like to add some subtle hillshading. I can see hillshading is built into some of the other styles like Satellite Afternoon.

How could I include the hillshade source to my OSM Bright 2 style?


texts: Performing row update within Select By Location based on field within selected feature class using ArcPy?

texts:
I am trying to select address points from address feature class within a selected zip code from zipcode feature class and populate the address feature class field "ZipCode" with the "Zip_Code" field from the zipcode feature class that is selected.  What I have so far doesn't produce an error, but continues to run and the ZipCode field is not being populated.

I am new to python.

Following is what I have so far.

#Set variables for feature classes
addresses = "Addresses"
zipcodes = "Zip_Codes"

#Perform Select by Location to create the zipcode data from ZipCodes to Addresses
#Create layer file
address_layer = "addr"
zipcode_layer = "zip"

if arcpy.Exists(zipcode_layer):
    arcpy.Delete_management(zipcode_layer)

arcpy.MakeFeatureLayer_management(zipcodes, zipcode_layer)
result = arcpy.GetCount_management(zipcode_layer)
print("There are {0} zip codes in the Zip_Codes feature class".format(result))

if arcpy.Exists(address_layer):
    arcpy.Delete_management(address_layer)

arcpy.MakeFeatureLayer_management(addresses, address_layer)

#Update zipcode field in Addresses feature class based on selectbylocation on Zip_Codes
#Initiate count to 1
count = 1
fieldListA = ["ZipCode"] #Zip code field from addresses variable
fieldListZ = ["Zip_Code"] #Zip code field from zipcodes variable
with arcpy.da.UpdateCursor(addresses, fieldListA) as upCursor:
    for row in upCursor:
        #Select all address points within the selected Zip_Code field within zipcode_layer
        arcpy.SelectLayerByLocation_management(address_layer, "WITHIN", zipcode_layer)
        result = arcpy.GetCount_management(address_layer)
        #print("There are {0} addresses".format(result))

        strZip = row[0]
        upCursor.updateRow(row)

        count = count + 1


texts:
I'm trying to get an elevation profile from a geotif file ( downloaded from http://srtm.csi.cgiar.org/). I downloaded the libgeotiff library but i can's see any function that can return an elevation for given pixel/coordinate or any useful documentation/example to how to do that.
Can someone help me out how to get the elevation? i want to do it in c\c++


texts:
I cannot get QGIS 3 to work. I've filed bug report #19629 but thought I'd also ask here to see if anyone has had this problem as it may just be my installation.

This is on Windows 10.
I have uninstalled all versions of QGIS and renamed the .qgis2 directory and then reinstalled v3.2.1 with the standalone installer.

Start QGIS 3.2.1
click the "Open Shapefile Layer" toolbar button
QGIS opens the "New Shapefile Layer" dialog
click the "..." button next to the File Name input
QGIS opens a "Save Layer As" dialog and immediately CRASHES to the "Oh Uh!" dialog

It also crashes if I start the program and click the "Open Project" button.
QGIS opens "Choose a QGIS Project File to Open" dialog and then crashes to the "Oh Uh!" dialog with no further input from me.

None of these crashes have anything to do with data. I never get to select a file.

QGIS v2.18 works but I'd like to try v3.   :)


texts: How to make button with html and mapfile?

texts:
I have several raster layers defined in mapfile:

   LAYER
     NAME "01"
     TYPE raster
     STATUS off
     OPACITY 30
     DATA "01.png"
  END

   LAYER
     NAME "02"
     TYPE raster
     STATUS off
     OPACITY 30
     DATA "02.png"
  END

   LAYER
     NAME "03"
     TYPE raster
     STATUS off
     OPACITY 30
     DATA "03.png"
  END

They are in .png format in raster folder on my server (data/www/project4/raster).

I need to insert button which will update my map with raster. First button is 01, and when I click on that button, raster should be rendered on my map. Also, I'm interested is there an option where I can click on some button just to open menu with three buttons which will lead rendering those rasters.

Thank you.


texts: How to sample raster datasets at ALL points?

texts:
I know it is possible to use tools like ArcGis's 'extract sample' tool to extract values at a given number of points. I would like to sample every pixel in a raster, and drill down through all other rasters (which have the same extent and resolution). How could I do this? I can use QGIS, Arc, or Idrisi.

Would the best way be to generate regular points at a spacing of 'RESOLUTION' and an initial offset of 'RESOLUTION/2'?


texts: How to show a geoTiff on an iOS application?

texts:
I'm developing an iOS application that will use offline maps. These maps are on GeoTiff format but I don't know how to do it.

I've found this question: How to display a GeoTIFF on mapView but it uses ARCGis and I don't have a license to use it.

Is there any other way to show a GeoTiff? And, of course, use GPS to show where user is.


texts: Lack of Options in Gdal Tool "Build Overviews" QGis

texts:
I am using QGis 1.8 and I want to build overviews (pyramids) for satellite imagery using the Gdal Tool. I would like to use the resampling method "cubic convolution", but the option doesn't exist although it is listed in the online help (http://www.gdal.org/gdaladdo.html). Two other options are in my list that don't show up in the online help however: "nearest" and "mode".

Does anybody know why "cubic" isn't there or can recommend me one of the other methods for use with satellite imagery?

Thanks a lot in advance,
Elena


texts: Calculate River Runoff Volume (Q) with GRASS

texts:
I want to calculate Hydropower between points at a River. I already created a Rivernetwork and the Accumulationmap (without weighting) with r.watershed.
My next Step was to set equivalent points along the River and find out their elevation -> done.
So now I want to calculate the Rivervolume for each point (along the River).
For this purpose I downloaded Precipitation data.

My Idea was to multiply the Accumulationmap with the Cell Size and the Precipitationmap (Previously, the card was processed that only the outflow remains on each cell (Runoff=Precipitation - Evapotranspiration)). But I´m not sure if that makes sense because the Accumulation says only how much Cells flow in but not which one (and there are Different Runoffs from different cells -> some also negative).

Do you have an idea how I can realize that and calculate the Riverrunoff for the points?


texts: Create cross sections equally spaced in a poly line imported from QGIS to MATLAB?

texts:
I had to create an Excel sheet with each profile cross-section which gave me the coordinates and elevations from a poly line of 2300 m length, referenced from 5 different DGM, which composed the study area. Each cross-section had different length and was spaced by eyelet roughly 100 m. I need now to space them by 5 m and I will get a really big data.

How can I create and export the cross sections profiles equally spaced and with the same length from the data above described?


texts: How do you make a function for changing projection with proj4js in Openlayers?

texts:
I have included Proj4js in OpenLayers. I would like to create a function (so I can use a make a button). To change the projection of my map.

I have included in javascript

Proj4js.defs["EPSG:54008"] = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";

And I have a added a GeoJSON layer.

        var world = new OpenLayers.Layer.Vector("World", {
            projection: new OpenLayers.Projection("EPSG:4326"),
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "data/worldcountries.geojson",
                format: new OpenLayers.Format.GeoJSON()
            })
        });

How can I make a function so I can change the projection of the vector-layer.

Hope someone can help


texts: Convert an arbitrary meta-data-free map image into QGIS project

texts:
My background: experienced programmer, amateur historian, complete beginner to GIS (so I'm using QGIS). I did learn to use the QGIS georeferencer tool though.

The input: suppose I have an image like this Wikipedia: File:VoyagesOfRabbanBarSauma.jpg. Assume that it's a "naked" JPG file without any geographic metadata and the author cannot provide any such information. The input image is licensed under GNU FDL, which I'm happy to release my derivative works under.

The goal: import the path shown in the Wikipedia map into QGIS so I can change the projection, superimpose elevation, add labels for physical geography, etc.

The caveat: one workaround might just be to enter the lat/long of the locations in the map into a new GIS project. I don't want to go this route because I'd like to have a record of that particular author's estimate of the geographic information.

What I've tried: I have applied the QGIS Georeferencer to the raster JPG file with a bunch of points around the continents' edges. The alignment of points in the interior of the continent are so-so.

What I'd like to know: the distortions of this particular map are pretty noticeable around the northern edge of Eurasia, and it seems like it's a pretty standard projection. It seems to me that there ought to be an algorithmic way to identify the projection and projection parameters that generated this map, and recover an accurate non-interpolation conversion between image pixels and lat/long?

I see two questions on this:

Identifying Coordinate System of Shapefile when Unknown?
Unknown Coordinate System on old drawing

These seem to talk about a brute-force method involving GDAL or QGIS, and to a GIS-beginner, seem somewhat labor-intensive (basically georeference and fit them to known projections) but I'm willing to try them and see if I can automate them.

The instructions on both these GIS.SE are a bit higher than my experience allows, and any hints on how exactly to do this brute-force method would be very appreciated. E.g., question 2 above shows how to use gdaltransform to convert a list of coordinates into lat/long. Are the input coordinates x-y pixel locations?


texts: How to transform data from a custom projection in JTS?

texts:
I have some data in NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet

I wish to programmatically transform the coordinates to Lat long in WGS84 datum using JTS.

I know that JTS can transform the coordinates using the following code:

import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;

CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:23032");

MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
Geometry targetGeometry = JTS.transform( sourceGeometry, transform);

But my projection does not have an EPSG code.

How do I define the coordinate system in this case? How do I transform my data into EPSG:3857?


texts: How to find the Shortest Path between two coordinates points using pgrouting A-Star?

texts:
In PostGIS I have two tables: an osm road network and a table that contains the start and end coordinate points of each path (i.e. userid, start_x1, start_y1, end_x2, end_y2). I would like to know:
1) how to create an osm road network that can be used for implementing the Astart algorithm with pgrouting (I actually have an osm road network which I obtained using osm2po and that works fine Dijkstra shortest path algorithm using pgrouting layer plugin)
2) how to find the shortest path of each start and end coordinate points using A-star.
Thank you, Selini


texts: For some reason one polygon in a shapefile is hidden behind another. How to solve this?

texts:
In this photo, you can see that Great Brittain consists of two features. London and the rest.

Problem is that London doesn't show as a separate colour. When I select great brittain, the area where London is supposed to be is also selected, but London in the attribute table is not. This shouldn't happen because I made this from a NUTS 1 map, where I gave London code 32 and the rest 31 and I dissolved on these numbers. And no, London is not in another layer.
However sometimes, when I zoom in and out quickly, when the rest of Great Brittain hasn't loaded yet, I see London:

Does anyone know how this happened and how to get London back, which seems to be behind GB for some reason?


texts: Why is my count in polygon number greater than the number of raw data points?

texts:
I'm working with GPS points and want to count how many of them fall within protected areas of 4 countries in Southern Africa (Eswatini, Zimbabwe, Mozambique and South Africa). I get the shape files for each from here https://www.protectedplanet.net/

The issue is there are 1768 GPS data points but I get a count of 1990 points within the protected areas which should be impossible.

I thought it would work out if I did it separately but no luck. The number of overlapping points are 59 from Eswatini, 819 from Mozambique,154 from Zimbabwe and 958 from South Africa which still sums to 1990.

Then I thought maybe some of the protected areas overlap across countries so some points are counted twice but it doesn't look like that's the case either.

I've carried it out in QGIS using Count points in polygon and in R with st_intersects and get the same result.

Here are the GPS data if you'd like to replicate what I've done https://drive.google.com/file/d/1IXTCQXAmQQehxpTmJn9aB4gfp1Bm2RH2/view?usp=sharing


texts: Showing OSM source acknowledgment on the map in QGIS?

texts:
I am using OpenStreetMap as a base layer via XYZ tiles in QGIS 3.01.

I can't seem to find a way of having source acknowledgement appear on the lower right corner of the map.

Can somebody point me in the right direction?


texts: Accentuated characters not displaying well in Arcpy Toolbox in ArcMap

texts:
Accentuated are not well displayed in my Toolbox made with ArcPy.

For instance :

def __init__(self):
    """Define the tool (tool name is the name of the class)."""
    self.label = "Sample éééé"

Produce the following output while executing from ArcMap :

Sample Ã©Ã©Ã©Ã©

é become Ã©


texts: What can I do with gdalwarp that I can't do with gdal_translate, and vice-versa?

texts:
I'm new to GDAL and its utilities.
My experience with both these programs (Warp and Translate) doesn't go much further than doing some image clipping, whether by coordinate's bounding box (that I do with both programs), "size"'s bounding box (that I only can do with translate) or use a polygon cut-line (that I only can do with gdalwarp).

Besides these clipping operations, what other operations can I do? And from those which are exclusive from each program.


texts: Why do I only get a blue background with TileMill and OpenStreetMap?

texts:
I use OpenStreetMap Carto in Tilemill in order to export mbtiles for a small town. I managed to setup everything according to the installation guide. However, when I open the OpenStreetMap Carto project in Tilemill I only get a blue background which is set by Map {background-color: @water-color;} in style.mss. Other than that there is no rendered map to be seen. I checked my database in QGIS and there it looks ok. Since I'm very new to this, I've no idea where to look for the problem. Any suggestions?

PS: I had a problem similar to this one: Error with TileMill and OpenStreetMap Carto project, and was able to resolve it by the suggested solutions the result being that the errors didn't show up but I ended up with the map showing only its blue background.


texts:
I'm working in the development of an user interface in which the user should be able to enter a UTM coordinate in the following format:

I don't know much about these type of systems, but I do want to provide a good user interface with client side validation. I know that the zone will be 1-60 with a plus or minus sign to determine the hemisphere. I would like to know what would be the valid range to input for easting and northing field. Currently I'm just checking if 'Double' number is provided, but I know the valid range should be more narrow.

Thanks